CONVEX ENTROPY, HOPF BIFURCATION, AND 
VISCOUS AND INVISCID SHOCK STABILITY 



BLAKE BARKER, HEINRICH FREISTUHLER, AND KEVIN ZUMBRUN 



Abstract. We consider by a combination of analytical and numerical techniques some basic ques- 
tions regarding the relations between inviscid and viscous stability and existence of a convex entropy. 
Specifically, for a system possessing a convex entropy, in particular for the equations of gas dynam- 
ics with a convex equation of state, we ask: (i) can inviscid instability occur? (ii) can there occur 
viscous instability not detected by inviscid theory? (iii) can there occur the — necessarily viscous — 
effect of Hopf bifurcation, or "galloping instability"? and, perhaps most important from a practical 
point of view, (iv) as shock amplitude is increased from the (stable) weak-amplitude limit, can there 
occur a first transition from viscous stability to instability that is not detected by inviscid theory? 
We show that (i) does occur for strictly hyperbolic, genuinely nonlinear gas dynamics with certain 
convex equations of state, while (ii) and (iii) do occur for an artifically constructed system with 
convex viscosity-compatible entropy. We do not know of an example for which (iv) occurs, leaving 
this as a key open question in viscous shock theory, related to the principal eigenvalue property 
of Sturm Liouville and related operators. In analogy with, and partly proceeding close to, the 
analysis of Smith on (non-)uniqueness of the Riemann problem, we obtain convenient criteria for 
shock (in)stability in the form of necessary and sufficient conditions on the equation of state. 
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1. Introduction 

Stability of shock waves has been the subject of intensive investigation over the more than 50 
years since the question was opened in the inviscid setting by Landau, Kontorovich, Dy'akov, Lax, 
Erpenbeck, and others in the early 1960's. Before that period, shocks seem to have been assumed 
to be universally stable, and in fact under most reasonable conditions, they are so. See |BE| for an 
interesting general discusion of these early investigations, and jS[ IMP] on the related question of 
uniqueness of Riemann solutions; we recommend also the original sources |ErH IEr2| IG] \M1\ \M2\ 
IM3j . Starting in the mid-1980's, these investigations have been widened to include also viscous 
shock stability, or stability in the presence of regularizing transport effects such as viscosity, heat 
conduction, magnetic resistivity, and species diffusion. See the surveys |Zll \Z2\ IZ31 IZ5j for general 
accounts of progress in this direction, and |BiBl IGMWZl] IGMWZ2] for accounts of progress on the 
related inviscid limit problem. See also the recent investigations [TZll ITZ2| ITZ31 ISSl IBeSZj on the 
possibility of Hopf bifurcation of viscous shocks, and stability of resulting time-periodic "galloping" 
shocks. 

At a technical level, stability and bifurcation of shock waves is now fairly well understood, both 
in one and multi-dimensions. In particular, nonlinear inviscid shock stability reduces [MlJ to a 
spectral condition that may be checked by evaluation of an associated Evans-Lopatinski function 
commputable by standard linear algebraic operations; see [Erl] IMll IT], IFP] for studies of the 
Lopatinski function in various contexts. Likewise, determination of nonlinear viscous shock stabil- 
ity/bifurcation has been reduced to verification of spectral conditions on the linearized operator 
about the wave, that may be readily checked by efficient and well-conditioned numerical Evans 
function computations as described in |Brt IBrZ} IHuZll \ZA\ IZ5] (see, for example, the recent studies 
[HLZl IBHRZl IBHZ11 iBLeZl iBLZl |HLyZl] |HLyZ2[ EH). (Cf. [mIZ3l iPZl IFST1 IFS2] . [HLZl Thm. 



1.4], |BHZH section 4] for situations in which the spectral conditions may be checked analytically.) 

However, at a level of basic intuition/understanding, some important questions remain, even in 
the one-dimensional setting. In particular, for all of the above-mentioned investigations of physical 
systems, whether inviscid or viscous, shocks were seen to be one-dimensionally stable. Indeed, 
as regards physically relevant systems of continuum mechanics, one-dimensionally unstable waves 
have so far been found only for media permitting phase transitions [GZ] \Z6\ IZMRS} IBM} [TZ4], 



and established wisdom [Erl[ IBE|, IMP| seems to say that instability is associated with effects not 
modeled in an ideal gas equation of state. 

One could ask, therefore, whether there might be some simple and commonly satisfied structural 
condition guaranteeing one-dimensional stability. In particular, in the simplest and most famil- 
iar setting of inviscid gas dynamics, could the classical condition of thermodynamic stability, or 
convexity of the equation of state 

e = e(r, S) 

relating internal energy e to specific volume r and entropy S by itself be sufficient to imply one- 
dimensional stability of shock waves ? Remarkably, the answer to this natural question up to now 
does not appear to have been knownj^J 

A natural generalization pointed out by Lax [L] of the property of thermodynamic stability to 
arbitrary systems of inviscid conservation laws 

(1.1) u t + f(u) x = 

is existence of a convex entropy, i.e., a function ij(u) satisfying 

d 2 r] > and drjdf = dq 

for an appropriate entropy flux function q(u), whence, for smooth solutions, 

V(u)t + q{u) x = 0. 



One could wonder whether for general systems (1.1), existence of a convex entropy is sufficient to 



imply one- dimensional stability of inviscid shock waves. This question includes the above one as a 
special case, as in the case of gas dynamics, one can take rj as the negative of the thermodynamical 
entropy S (r, e) , the Legendre transform of e 
For viscous systems of conservation laws 

(1.2) u t + f(u) x = (B(u)u x ) x , 

shock waves are heteroclinic traveling waves. It is known |GZl IZS| that viscous stability is a 
stronger condition than inviscid stability. However, up to now it is not known whether this logical 
implication of inviscid by viscous stability is strict for any interesting general class of constituents 
(/, B). Destabilization by viscous effects of an inviscidly stable shock wave would be a physically 
interesting phenomenon if it occurs, and similar phenomena arising in Orr-Sommerfeld theory and 
incompressible flow make this not implausible. On the other hand, if this could be shown not to 
occur, that would be equally interesting, and would greatly simplify the verification of stability, 
reducing this to the study of the simpler inviscid Lopatinski function, a linear-algebraic quantity, 
rather than the viscous Evans function, defined in terms of solutions of an associated eigenvalue 
ODE. In this connection, we were wondering whether the existence of a viscosity-compatible entropy, 
n as above with now also 

d 2 rj B>0, 

might imply equivalence of viscous and inviscid instability. 

As pointed out in |Z1[ \Z2\ \TZ2\ ITZ3j , the same analysis [GZ1 IZSj showing that viscous stability 
implies inviscid stability, shows that, as long as the viscous profile remains transverse as the inter- 
section of the invariant manifolds as which it is defined, viscous destabilization of a stable inviscid 
shock necessarily must occur through the passage through the imaginary axis into the unstable 
(positive real part) half plane of a complex conjugate pair of eigenvalues of the linearized opera- 
tor about the wave, a leading-mode Hopf bifurcation. Here, we are imagining a transition from 
viscous stability to instability as some bifurcation parameter, typically shock amplitude, is varied. 



^The suggested counterexample of [G] is incorrect, as we show below; see Remark 
2 5* is concave if and only if e is convex. See [MP] or Lemma 6.7 in [Zl] , 
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"Galloping instabilites" arising through leading-mode Hopf bifurcations are familiar in detonation 
theory [TZ4] , but have up to now not been observed in the shock wave context. Could it be that the 
existence of a viscosity-compatible entropy precludes complex, or even just purely imaginary eigen- 
values? Or could it possibly imply a "[strong] principal eigenvalue property" that any non-stable 
eigenvalue, of the linearized operator about the wave, with largest possible real part must be real 
[and simple]? In both cases, leading-mode Hopf bifurcation would be impossible. Could it even be 
that existence of a viscosity-compatible entropy implied both transversality of the shock profile and 
impossiblity of leading-mode Hopf bifurcation? In that case, viscous and inviscid stability would 
coincide. 

In this paper, we examine these and related questions using a combination of analytical and 
numerical techniques. The following are the main results of this examination: 

Theorem A (on gas dynamics), (i) There exist equations of state e = E(t,S),t > 0, S € K, 



that (a) satisfy the standard assumptions (J1)-(J4), (G1)-(G6), and (H1)-(H4) detailed below, with 



in particular, the largest and smallest eigenvalues simple and genuinely nonlinear, entropy and shock 
speed monotone along the forward 1- and 3-Hugoniot curves, and a global concave entropy function 
S(r,e) defined onr,e> 0, and (b) admit inviscidly unstable shock waves, 
(ii) The equation of state 

(1.3) e -(r,S) = -+CV/ c2 -/ c , C » 1, 

r 

is an example for (i). 



Numerical Observation A (on gas dynamics). For the equation of state (1.3), in all cases 
we investigated numerically, (a) viscous [in] stability is equivalent to inviscid [in] stability, 

(b) the viscous- stability problem has no non-zero imaginary eigenvalues; in particular, transitions 
from stability to instability occur exclusively by real eigenvalues passing through the origin, and 

(c) in situations of instability, the eigenvalue with largest real part is real and simple. 



Numerical Observation B (on general systems). There exist viscous systems (1.2) of 
conservation laws, endowed with a compatible entropy, that 

(a) admit shocks that are inviscidly stable, but viscously unstable, 

(b) the viscous- stability problem sometimes does have non-zero imaginary eigenvalues, while 

(c) in all situations of instability we investigated numerically, the eigenvalue with largest real part is 
real, and transitions from stability to instability occur exclusively by real eigenvalues passing through 
the origin, 

(d) in some cases, there are an even number of unstable (and all real) eigenvalues, and 

(e) in some cases the eigenvalue with largest real part is not simple. 

Our analysis divides roughly into an inviscid and a viscous part, which we now sketch. 

1.1. Inviscid analysis. In the first part of the paper, we revisit the inviscid stability analysis 
of Erpenbeck-Majda |Erl| |MT| |M2"1 IM3] and reconcile it with the analysis of Smith 0, IMP] on 
nonuniqueness of Riemann solutions, a related but slightly weaker condition than instability of an 
individual component shock. Let 

P = P(t, e) := -e T (r, S(r, e)) 

denote the pressure function determined by inversion of the equation of state e = e(r, S) for r 
fixedEl 



pressure law p = p(r, e) may also be considered in the absence of an equation of state, being sufficient by itself 
to close the equations of gas dynamics |MP ( Sm . For the observation that any pressure law can to some (interesting!) 
extent be interpreted as being associated with an equation of state, cf. Sec. 3.3. 
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Smith formulates criteria on the equation of state, amounting to a weak condition 

(Weak) — _~ T _ S < — — , or, equivalently, p T < 

e s e TT e T 2 

a medium condition, 

_2 

e e T _|_ ^ 2 
(Mediumu) — zz^z — < — 2eeTT _ , or, equivalently, p T < — , 

and a strong condition, 

e TS 1 

(Strong) — -z-^z. — < — — , or, equivalently, p T < 0, 



and derives in particular that (Mediumu ) is equivalent under mild additional hypothese^jto unique- 
ness of Riemann solutions. 

Recall now some standard properties of a pressure law p = p(r, e), r, e G M + : 



(Jl) 
(J2) 
(J3) 



V — v( T i e ) > 0- (Positivity) 
(d T — pd e )p < 0. (Hyperbolicity) 
(d T — pd e ) 2 p > 0. (Genuine nonlinearity) 



(J4) p e > 0. (Weyl condition) 

The following is the technical centerpiece of Part I of this paper. 



Theorem 1.1. Under hypotheses (G1)-(G6), (|Hl|) -(|H2| ) of Smith (see Sec. 3.2 below) on the 
equation of state e, or, more generally, assumptions (J1)-(J4) on the pressure law p, positivity of 
the signed Lopatinski determinant (see Sec. 2.1) (for all shocks) is equivalent to 



(Mediums) 



< 



yJ1ee TT 



+ 1 



or, equivalently, p T < cp/v2e = p 



PPe - Pr 



2e 



while ( |Sm| :) uniqueness of Riemann solutions (for any data) is equivalent to (Mediumu). The 
four conditions are related by 



(Strong) (Mediumu) (Mediums) =>• (Weak). 



In particular, condition (Strong) by itself is sufficient to imply stability of all shocks, while violation 
of (Weak) implies existence of unstable ones. 

Theorem A is a corollary of Theorem 1.1 and the following finding. 



Theorem 1.2. Equation of state (1.3) satisfies (G1)-(G6), (H1)-(H2) ; and violates (Weak). 



Both the particular implication ( Mediumu ) => ( Mediums ) in Theorem 1 . 1 and the general mean- 
ing of the signed Lopatinski determinant are elucidated by the following fact, which holds for 
arbitrary systems. 

Theorem 1.3. If the signed Lopatinski determinant A(a) of a family (U-(a), U+(a)),a 6 M. of Lax 
shocks undergoes a sign change at some critcal value a* G IR, the Riemann problem loses uniqueness 
near the initial data (Ui,U r ) = (U-(a*),U+(a*)). 



4 Cf. Theorem 1.1, within which we include Smith's result for the reader's ease. 
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Remark 1.4. A general polytropic equation of state e(r, 5) = e s / Cv /r r , r,c^ > 0, satisfies (Strong) 



or T 2 < T(r + l). By Theorem [D] we see that a local counterexample e(r,S) = e s /r + f(S), 
f » 1, proposed by C. Gardner (final paragraph of [G]) is incorrect, as this equation of state 
likewise satisfies ((Strong), by essentially the same computation. Nonetheless, Gardner's larger 



assertion, that "... all of the situations we have considered are logically possible, without violating 
the fundamental condition for local thermodynamic stability, namely that e is a convex function of 



r and S" turns out to be correct, as shown by example (1.3). 

Remark 1.5. Taylor expanding ( |1.3[ ) in C, we obtain also a simplified local example 
(1.4) e(r,5) = e s /t + S - Ct + t 2 /2, C » 1, 

of a convex equation of state permitting instabilities, for which Hugoniot curves, monotonicity of 
S and ex, etc. may be computed explicitly, giving a concrete illustration of the theory. This does 
not satisfy 



Proposition 



Jl), but may be treated by similar techniques as used to establish Theorem A; see 



4.1 below. On the other hand, we find that the seemingly similar example 
e(r,5) = e s /T-Cr + T 2 /2, C » 1, 



does not permit instabilities, despite violating (Weak); see Proposition 4.3 



Remark 1.6. Example (1.3) is the more surprising in view of recent results [LVj in nearby set- 
tings showing that existence of a convex entropy implies stability. Indeed, Theorem A appears 
at first sight to contradict Theorem 2 of [LVj . which asserts that, under some mild technical as- 
sumptions, existence of a convex entropy, simplicity of extremal characteristics, and monotonicity 
of entropy and shock speed along forward Hugoniot curves together imply inviscid stability of 
arbitrary-amplitude extremal (i.e., 1- or 3-) shocks. This would be very interesting to resolve; at 
the least, it shows that there is a very narrow window between negative and positive results. 

Remark 1.7. Smith's weak and strong conditions are phrased in terms of the gas law e = e(r,p) 



obtained by inverting p = p(r, e) with respect to e, under the assumption that p e > (( J4) below). 
Our versions are equivalent to those of Smith if and only if p e > 0, through the relation e T = 
but remain valid also in the general case, as Smith's therefore do not. Similar observations are made 
in |MPj regarding the logical ordering of Smith's conditions, which again requires p e > 0. Thus, our 
conditions above are in fact extensions of Smith's conditions into the realm p e < 0, or, equivalently 
(by ps = p e {de/dS) = p e T), ps < 0. In this case, other aspects of Smith's global theory break 



down, in particular, the conclusions of Proposition 3.2 and, thus, Theor em |1.1| However, we find 



that individual 1-shocks ([/_,[/+) are stable for p e \u+ < 0; see Corollary 2.6 



1.2. Viscous analysis. In the remainder of the paper, we investigate viscous stability of the 
above and other systems via numerical Evans function computations, seeking viscous instabilities 
not predicted by inviscid theory, and in particular purely imaginary eigenvalues of the linearized 
operator about the wave. The Evans function, as described for example in |AGJ|. IFW] IGZ| . is an 
analytic function D(\) associated with a viscous shock wave, defined on the nonstable complex 
half-plane {A : 3ftA > 0}, whose zeros encode the stability properties of the wave, in particular, 
corresponding away from the origin with eigenvalues of the linearized operator about the wave. See 
\Z1\ IZ21 IZ31 IZ51 ITZ21 ITZ3} ISS} IBeSZ] and references therein for more detailed discussions of the 
relation between the Evans function, spectral, linearized, and nonlinear shock stability, and Hopf 
bifurcation. 

As described in |BrZ} IHuZl} IZ4| IZ5| , the Evans function may be efficiently computed by shooting 
methods as a Wronskian evaluated at x = of decaying modes at x = ±oo of the eigenvalue 
ordinary differential equation (ODE) associated with the linearized operator about the wave, and 
this has by now been carried out successfully for a number of interesting systems/ waves arising in 
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continuum mechanics using the standard Runga Kutta 4-5 (ode45) routine supported in MATLAB 
in conjunction with the exterior product or polar coordinate algorithms |Brt IBrZl IHuZlj of the 
MATLAB-based Evans function package STABLAB [BHZ2] , in particular for the computationally 
intensive problem of stability of detonation waves |HuZ2l IBZH IBZ2} |BHLyZ2|). 



However, the equations of state (1.3) and (1.4), due to the separation of scales introduced by 



the large parameter C » 1 and exponential dependence on parameters, induce a computational 
difficulty far beyond any of the systems so far considered. In particular, computations with RK45 
were not practically feasible, even on a parallel supercomputing machine. To carry out Evans 
function computations for this system, rather, we found it necessary for the first time to use 
an ODE solver designed for stiff systems, namely, the odel5s routine supported in MATLAB, a 
recently-developed algorithm based on modern numerical differentiation formula (NDF) methods. 



To indicate the relative stiffness of the systems associated with (1.3)— (1.4) as compared to that 
of previously considered cases, the performance of Evans computations with ode45 and odel5s are 
similar for all previously computed examples; however, in the present case, odel5s outperforms 
ode45 by 2-3 orders of magnitude, yielding performance that is not only feasible, but in the general 
range seen for nonstiff computations. 

The results of these numerical computations, gathered above as "Numerical Observation A (on 



gas dynamics)" hold also for equation of state (1.4). 



The other results, presented above as "Numerical Observation B (on general systems)", are 
obtained by computations on an artificially constructed 3x3 system of viscous conservation laws 
(1.2 ^jof a form suggested by Stefano Bianchini [Bi] , Specifically, as shock amplitude is increased, 



several eigenvalues cross the origin into the unstable half-plane, some of them eventually coalescing, 
splitting off of the real complex conjugate pair, and crossing back into the stable half- 

plane through the imaginary axis at other points than 0. This is a Hopf bifurcation, though not 
a leading- mode one. Note finally that the examples with an even number of unstable eigenvalues 
exemplify a scenario not detected even by signed versions of the inviscid stability condition, — i.e., 
essentially, the signed Lopatinski determinant — that find the parity of the number of unstable 
roots as described in |GZ| IZ1 [ IZ2j. 

1.3. Discussion and open problems. In conclusion, we have confirmed the claim of C. Gardner 
[G] made almost 50 years ago, but since then apparently not paid much attention to, that the 
equations of gas dynamics, even under all of the usual thermodynamical and structural assumptions 
associated with a "normal" gas, can support unstable inviscid shock waves. Moreover, we have shown 
this in a somewhat stronger sense than that apparently envisioned by C. Gardner, demonstrating 
that this can be accomplished not only locally, but globally, with a convex equation of state e = 
e(r,S) satisfying expected asymptotics as S — > ±oo or r — > 0,+oo. In the process, we shed new 
light on the important work of R. Smith [S] and Menikoff-Plohr |MPj on uniqueness of Riemann 
solutions. 

This in a sense validates the inviscid stability theory of Erpenbeck-Majda, showing that stability 
does not hold "automatically" for equations of state satisfying all conditions of a classical gas. On 
the other hand, we find so far no distinction between viscous and inviscid stability for shock waves 
in gas dynamics. We do, however, show by example that such a distinction, and in particular the 
new phenomenon of Hopf bifurcaton, can occur for a 3 x 3 system with global convex entropy. 
Whether such phenomena arise also for Lax shocks of physical systems such as gas dynamics, 
magnetohydrodynamics (MHD), or viscoelasticity, must be left to further investigation. We point 
in particular to the viscous-stability analysis of full (nonisentropic) MHD, so far not carried out 
even for the classical polytropic equation of state, as a promising candidate and an important 
problem for further study. 



The same size as the equations of gas dynamics. 
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Regarding the initial transition from stability to instability, we find no distinction between viscous 
and inviscid behavior for any of the systems we consider. Likewise, we find no counterexample to 
the (weak) "principal eigenvalue property" that the largest unstable eigenvalue of the linearized 
operator about the wave be real. It would seem most interesting to know whether this property 
holds in any generality. 

Finally, recalling that shock instability in the physical literature has often been associated with 
the phenomenon of phase transition |Erl4 IBE[ IMP] , we repeat that our examples here do not 
involve phase transition; rather, the principal mechanism for instability, at a technical level, i.e., 



the feature leading to violation of (|Weak|), is stiffness, or appearance of multiple scales via the 



parameter C » 1. Indeed, computing p(r, S) = —e T (T, S) = ^ + C — r for example (1.4) reveals 
a similarity to a "stiffened" or "prestressed" equation of state sometimes used to model shock 
behavior near the liquid state. For example, water is often described by 

(1.5) p = Tpe - 7P0 

with r and the base stress Po determined empirically (ITj, IHVPM) ; similar techniques are used to 
model liquid argon, nickel, mercury, etc. fHl lCDM ]. The careful numerical study and categorizaion 
of behavior of these and other more exotic equations of state such as van der Waals, Redlich-Kwong, 
etc., both with and without viscosity, appears to be another important direction for further study, 
and one in which surprisingly little has so far been done. In particular, a viscous counterpart of 
the inviscid analyses of stability carried out for van der Waals gas dynamics in [ZMRSj IBM) would 
be a valuable addition to the literature. 

PART I. INVISCID STABILITY FOR GAS DYNAMICS 

Sections 2 and 3, the first two of the three sections of this first principal part of the paper are 
devoted to establishing a network of algebraic and geometric conditions on both individual shock 
waves and Hugoniot curves that alone or together characterize (un)stable shock waves or imply 
or preclude the existence of such waves in the sense of necessary or sufficient conditions. Part I 
culminates in Section 4 that shows the existence of inviscidly unstable shock waves in gas dynamics 
with certain equations of state. 

2. Signed Lopatinski criterion and relations to other conditions 

In this section, we define the signed Lopatinski criterion, evaluate it for gas dynamics with very 
general equations of state, and characterize its connections with other conditions on the equation 
of state. 

2.1. The signed Lopatinski determinant. We begin by recalling the Lopatinski condition of 
Erpenbeck-Majda |Erl| IMH IM2} IM3j for stability of inviscid shock waves. Consider a general 
system of conservation laws 

(2.1) f°(w) t + f(w) x = 0, iiiGf 

and a Lax p-shock |Sm[ ISe2|, ISe3j U± with speed a, satisfying the Rankine-Hugoniot conditions 

(RH) a[f] = [/], 

hyperbolicity 

(2.2) a(A ± )real, A ± := (f^)' 1 f w (U ± ), 
of endstates U± and the Lax conditions 

a7 < • • • < a~ 1 < a < a~ < ■ ■ ■ < a~, 

(2.3) 1 - - p-l p - n 

a{ < ■ ■ ■ < al < a < aJ. x <•••<<, 



where a- are the (ordered) eigenvalues of A± := (/°) 1 f w (U- 



± 



In this setting, one-dimensional linearized and (under additional mild structural assumptions) 
nonlinear inviscid stability is addressed (see |M1[ IMZTl [MZ2, BSJ) by the Lopatinski condition 

(2.4) 5 := 6(U-,U+) := det(^° rf, . . [f% A° + r+ +1 , . . .,A° + r+) ? 0, 

where rj, j = 1, . . . ,p — 1 are the "outgoing" modes of A_ := (fz.) fvi(US) relative to the shock, 
i. e., a basis of the invariant space associated with eigenvalues aj such that aj — a < and r~T, 
j = p + 1, . . . ,n are the "outgoing" modes of A + := (/^,) _1 f w (U + ), relative to the shock, i. e., a 
basis of the invariant space associated with eigenvalues such that — a > 0, and A± := f®. 
Implicit in this description is the defining property of a Lax p-shock that the outgoing eigenvectors 
on the left (right) indeed number p — 1 (n — p). 

Remark 2.1. Note that we have allowed here by the inclusion of /° a somewhat more general 
formulation than in (1.1). In the simplest case f°(w) = w occurring in many applications, the 
computations considerably simplify. However, as we shall see in the full gas case, there can in some 
cases be an advantage in being able to choose coordinates more flexibly. A short proof of the case 
f(w) = w may be found in |ZSj . Redefining w := f°, we obtain A + = fw = /m(/m) _1 5 yielding the 
result by the observation that eigenvectors of A + are equal to A^r^. See also |Se2t ISe31 IMlj . 



Definition 2.2. Consider a given p-shock (U-,U+) for a system (2.1). Assume that there is a 
continuous one-parameter family (U~(a), U + (a)), < a < 1, of p-shocks with U^(l) = U± and 
U ± (0+) = U* for some state £/*, that the bases {r^~, . . . , {r^ +1 , . . . , r+} and thus 

5{a) := 5{U~(a),U + {a)) 
are (chosen) continuous in a, and that the expression 

(2.5) A = 5(l)/limsgn<5(a)), 

makes sense and has the same sign for all such homotopies. Then A = A(C/_, U+) is called a signed 
Lopatinski determinant of (U-,U+). 

Obviously, the signed Lopatinski determinant A is well defined ( — more precisely, while the 
Lopatinski determinant 5 is defined up to a non- vanishing factor, the signed Lopatinski determinant 
A is defined up to a positive factor — ) at least whenever 

the system's state space is simply connected, 

(2.6) a p is an everywhere simple, genuinely nonlinear eigenvalue, and 

all p-shocks with arbitrary fixed left state U- group as regular directed curves S P (U-). 

In the current Part I of this paper we use the signed Lopatinski determinant to detect instability 
according to the following evident fact. 



Corollary 2.3. If a given hyperbolic system (2.1) satisfying (2.6) admits a p-shock with A < 0, 
then it also admits a p-shock with 5 = 0. 

Proof. For a corresponding homotopy, obviously 

A(U~(a), U + (a)) > for sufficiently small a > 0, 

and the existence of an a* with A(U~ (a*), [7 + (a*)) = follows by continuity. □ 

For the (deeper) meaning of the signed Lopatinski determinant in the viscous context, see Part 
II, Sec. 6. 
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Remark 2.4. As noted in |Fo] in the context of gas dynamics, and |Sel| more generally, a change 
in the sign of A may also be interpreted in terms of multi-dimensional inviscid theory, signalling a 
transition from weak multi-dimensional stability in the sense of Majda (Mil IMH lM3] to exponential 
instability; see also |MP| \Z1\ IZ3| . 

2.2. The Lopatinski condition in gas dynamics. Specializing to gas dynamics, we now com- 
pute the Lopatinski determinant and signed Lopatinski determinant explicitly. 

Isentropic gas dynamics. The isentropic Euler equations in Lagrangian coordinates are 

Tt — v x = 0, 

(2-7) 

y ' v t + Px = 0, 

with a given pressure law p = p(r), where r denotes specific volume, v velocity, and p pressure. 

From A± = ^ , we readily find that = —c,c, c := \f—p T (j), and r f = a ±^j , so 

that hyperbolicity corresponds to p' < 0. 

Moreover, (RH) gives <t[t] = —[v], a[v] = [p], yielding a 1 = — [p]/[r], whence, for a 1-shock, 

a = -\/-H7FL and so we have t/°] = M ( ^ -[pi/ir] ) ' Thus ' 

5 = [t] det f 1 M = [r](<r - c) < and A > 



for all Lax shocks. Notice that condition (2.6) is satisfied trivially. Note that [r] 7^ 0, since otherwise 
[v] = by (RH), and there is no shock. 

Full gas dynamics. In Lagrangian coordinates, the full (nonisentropic) Euler equations are 

Tt ~ V x = 0, 

(2.8) v t +p x = 0, 

{e + v 2 /2) t + (vp) x = 0, 

with a given pressure law p = p(r, e), where r denotes specific volume, v velocity, e specific internal 
energy, and p pressure. We assume here that the pressure law originates from a thermodynamic 
equation of state 

(2.9) e = e(r,5), 

where S is entropy, through the relation p = —e T , assuming that temperature T = e~s is positive, 
so that (2.9) may be solved for S = S(r, e) to obtain p(r, e) = e(r, S(r, e)). 
Alternatively, for smooth solutions, we have the simple entropy form [SmJ 

n - v x = 0, 

(2.10) v t + Px = 0, 

s t = o,. 

where p = p(r, 5) = — e r (r, S) for equation of state (2.9). Choosing w = (r, v, S) T , we thus have 





A± 




and, essentially by inspection, 



(2.11) a\ = -c, a 2 = 0, a 3 = c, c := y/—p T , 
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, so that hyperbolicity corresponds to p T 



, we obtain 



< 0. 




Next, from (RH), we have, taking without loss of generality v 
[v] = -a[r], [e + v 2 /2]-- 



0, 



(2.12) 



so that, for a 1-shock (recalling the relation p T 





or 



(2.13) 6 = [T]( Ps c[p}+Tc 2 (a-c)), 

where all quantities are evaluated at [/+, and c := \J—p T as above denotes sound speed. Note that 
[t] / 0, or else [v] = 0, [p] = 0, and thus [e] = 0, all by (BE). But, [e] = [r] = and [r] = by 
es > implies that [S] = and there is no shock. 

Note that, in the decoupled case ps = 0, we get 5 = [t](Tc 2 (ct — c)) < 0, by T > 0, essentially 
reducing to the isentropic case. 

Note 



We assume, here and henceforth, that condition (2.6) holds for the Euler equations (2. 



that whenever Smith's conditions (G1)-(G6), (H1)-(H2) (cf. Sec. 3.2) hold, this is automatic. 



Theorem 2.5. For equations (2.8) with p T < a Lax 1-shock satisfies the signed Lopatinski 
condition A > if and only if 

(Lop) p s \p] < Tp T (a/c - 1), 

or, alternatively, 

(Lop a i t ) p T [p] < -co, 

where all quantities not in [.J are evaluated at U+. 



Proof. By homotopy from the general case to the decoupled case. Equivalence of (Lop a it|) and 
-c 2 =p T -pp e and pe(r,e) = ps(r, S)/e s (r, S). 



( Lop ) follows by p T 



□ 



We state the following elementary result mainly in order to illustrate the competition between 
isentropic and nonisentropic effects. 

Corollary 2.6. For a gas with equation of state satisfying T, p > 0, p T < 0, and ("anti-Weyl" 
condition) ps\u+ < 0, all Lax shocks satisfy the signed Lopatinski condition A > 0. 



Proof. In this case, each of the terms psc[p] and —Tp T (a — c) on the righthand side of (2.13) are 
strictly negative, by T > and p T < (hyperbolicity), so that sgn 5 = sgn[r]. But [r] < for any 
1-shock and [r] > for any 3-shock. □ 

Remark 2.7. Though we have restricted for clarity to the case that the pressure relation derives 
from a complete equation of state e = e(r, S), we could equally well carry out the analysis within 
the framework (2.8), with p = p(r, e), to obtain the condition p e \p] < c 2 (l — <r/c) or, more simply, 
(Lop a i t ); see Appendix [Aj This could alternatively be obtained indirectly, by the observation that, 

n 



for any choice of positive temperature function T = T(r, e) denned in the vicinity of U+, we can 
invert the Thermodynamic relation de = TdS — pdr to obtain a local entropy S and equation of 
state e = e(r, S); see Section 3.3 



Remark 2.8. Though the pressure law p = p(r,e) is sufficient to close the Euler equations, the 
Navier-Stokes equations, involving the Fourier law, require also the specification of a temperature 
law T = T(t, e). Both can be obtained from an equation of state e = e(r, 5), as described above. 



Relation to Majda's condition. We now verify that, assuming the Weyl condition p$ > 0, 
our condition (Lop) agrees with Majda's ID stability condition (more precisely, the signed version 
obtained by taking into account multi-dimensiional considerations; see Remark B.6, p. 523-524, 
Appendix C, [Zl]): 



(2.14) 



1 + M 



> 



T 



where T is the Gruneisen coefficient T := ^r-, and M 2 



1)M 2 , 



— , the subscript E denoting Eulerian 

C E 



values. Evidently, T > if and only if ps > 0. 

Rewriting (2.14) as jqj-^ < _^ M 2 using positivity of T and substituting ^£ for T, we have 

-[t]ps < 1 + M 



M 2 



We may compare to our condition (Lop) by expressing the (Eulerian) Mach number M in terms of 
Lagrangian quantities. Evidently, we must verify the identity — — — — — — — 



c[p] 



or 



1 + M c(c-a) 



M 2 a 2 ' 

In Eulerian coordinates, we have cte[p\ = [p v ]> or °~E = jf^j = = c, in agreement with the 
Lagrangian value. Moreover, v = [v] = — a[r] gives v — a = — c([t] + 1), so that, in case tq = 1, 
we have v — a = —o~t. On the other hand, ce = tc, hence M = and the desired identity is 



l + (v—a)/rc 
{v-<j)2/{tc) 2 



c(c—a) 



, or 



rc+(v—a) (c— it) 

(v-a) 2 



Without loss of generality (by rescaling) taking r_ 



1, 



-err, this becomes 



^ c 2 - , which is evidently true. 



Remark 2.9. Similarly as in Remark 1.7 we find that our condition (Lop) is an extension of Majda's 
condition into the realm ps\u+ ^0 where Majda's condition no longer applies. Indeed, for this case 
r < and (2.14) fails always, whereas we see by (Lop) that stability in fact always holds. 

2.3. Sufficient global conditions for (in)stability. The following sufficient conditions for sat- 
isfaction and violation, respectively, of the signed Lopatinski condition are as easy to get as useful. 
For definiteness, we restrict to the case of 1-shocks. 

Proposition 2.10. For a Lax 1-shock with [r] < 0, the signed Lopatinski condition A > is 
implied by 

1 



(Strong' 



The opposite, A < 0, is implied by failure of 



\u=u+ < 



IpY 



(Weak') 



ese TT U U+ < [p] ' 
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Proof. Because of [r] < 0, [p] is positive, and dividing both sides of (Lop) by the positive quantity 
c[p], we obtain 



( L °Ply 



< 



■f + 1 



+ 1 



[p] 



Recalling, (2.11), that the characteristic speeds for (2.8)-(2.10) are — c, 0, c, we have for a Lax 
1-shock that at the right state U + , — c < a < < c, or a < and |<r| < c. Hence, 1 < — ^ + 1 < 2, 
from which we obtain with (Lop 1 ), the stated stability conditions. □ 



Remark 2.11. Conditions (Strong') and (Weak') are global versions of conditions (Weak) and 



( Strong[) of [S]. Evidently, (| Weak[) =4> (|Weak'[) and (|Strong| =4> (|Strong~ 



2.4. Monotonicity of the Hugoniot curve. Using the formulation (Lop), we show how stabil 



ity of shocks relates to geometric properties of the Hugoniot curve, recovering both the classical 
observations of Erpenbeck and Gardner [Erl] [G] relating stability to geometric properties of the 
Hugoniot curve and the conditions of Smith for uniqueness of Riemann solutions. It remains to 
investigate monotonicity of the 1- (or, equivalently, the 3-) Hugoniot curve. From the assumption 
e~s = T > 0, we may solve e = e(r, S) to obtain S = ,§(r, e), and thereby 

(2.15) p = p{r,e) := -e r (r, s(r, e)). 
Combining the three equations ( RH| , we readily obtain the single equation 

(2.16) H(r,e) = [e] + (l/2)(p + p_)[T] = 



for (e, r) = (e+, r+), implicitly determining, together with ( |2.15 ), e as a function of r or vice versa. 
Alternatively, using (2.12) to relate different variables, we could consider this as determining p as 



a function of v, or r as a function of s, whose derivatives in each case may be computed by a 
straightforward application of the Implicit Function Theorem. We record the results here. 



Lemma 2.12. We have the change of variables formulae 



(2.17) 



1 



e s ' 



Pt 



+ 



Pe 



e s ' 



and, along the Hugoniot curve H = 0, viewing e and p as functions e = ej^(r) and p = Ph(t), 



(2.18) 



de H -\{Pt[t] +p + p- 



dr 



dpH 
dr 



Pt+P 



„ de H 



1 + ^ ™ dT 1 + ipeM 

Alternatively, considering v as a function of p along H = 0, we have, taking V- = 0, 



(2.19) 



2v 



dv 
dp 



[r] 



a 



P]Pe 



Proof. Relations (|2.17|) and (|2.18|) follow directly from the Implicit Function Theorem applied to 

f, den _ -|(p+P-)p e +Pr 



(2.15) and (2.16), where the final equality in (2.18) follows from — H 



dr 



together with c 2 = — p T + pp e . From v = — \J — [p] [t] (a consequence of (RH) and -u_ = 0, 
through [v] = — o~[t]), we obtain, differentiating implicitly v 2 = [p][r], — 2t>^ = [r] + [p] dp ^/ dr = 



^dpu/dr =. verifying by ( [238 ) 



□ 



Carrying out these computations, one may obtain various necessary or sufficient geometric con- 
ditions on the Hugoniot curve. In particular, we have the following result. 
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Corollary 2.13. Monotonicity, dv/dp < at a point [/+ on the 1-Hugoniot curve through U- is 
equivalent to 

„2 



(Monotone) 



e rS < 



1 



ese rr [p] 
Moreover, we have the string of implications 
(2.20) 



or, equivalently, p e \p] < cr 2 + c 2 . 



( Strong' ) ( Monotone ) ( Lop 1 ) => ( Weak' ) . 



Proof. From (2.19), we have evidently (by v, [r] < 0) that dv/dp < is equivalent to 

/dr — a 2 



> 



dp H /dr 



2 + I 



or, rearranging, p e [p] < a +tr as claimed, yielding (Monotone) by [p] > 0. The logical implications 
(2.20) then follow by the Lax condition \a\ < c. □ 



Remark 2.14. The implication (Monotone)=^ (Lop! ) is due to C. Gardner [G) . A slightly less sharp 
condition pointed out by Erpenbeck |Erl] is t hat m onotonicity of r with respect to s along the 
forward Hugoniot curve, ps < (see Remark 3.4), or 



(2.21) 



2<r 2 

e rS -gr 



e 5 e 



also implies (Lop 1 ), by \<r\/c < 1; likewise, (2.21) implies ([Monotone) and uniqueness of Riemann 
solutions. For other geometric implications of the various conditions, see Theorem 4.5 |MP| . 



Remark 2.15. From (2.20), we obtain ( Strong ) ( Lop-L ), and ( Strong )=>( Strong' ). 



3. Local conditions for (in) stability and the assumptions of Smith 

The above-obtained global conditions require knowledge of the Hugoniot curve in order to com- 
pute, which is in practice quite limiting. Adapting the ingenious observations of Smith, we now 
show that, under some simple structural assumptions on the Hugoniot curve, these global condi- 
tions may be replaced by local conditions more convenient for analysis. In this section, we establish 
Theorems 1.1 and 1.3. While the latter is proved in Subsection 3.4, Theorem 1.1 is an immediate 
consequence of Propositions 3.1 and 3.2 and Corollary 3.6. 

3.1. Local conditions and their requirements. We identify the following properties, satisfied 
for most standard equations of state. Denoting as the backward 1-Hugoniot through U + , H[(U + ), 
the set of all left states U- connected to by a Lax 1-shock, and assuming that H[(U+) is a 
directed curve, we distinguish 

(PI) [r] <0 on H[(U + ), 



(P2) p — > as U progresses along H[(U + ), 

(P3) e — > as U progresses along H[(U + ), 

(P4) r is increasing and p decreasing along H[(U+). 



Recalling now conditions (Strong), (Mediumyj), (Mediums), (Weak) of the introduction, we have 
the following detailed local characterizations of stability. 
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Proposition 3.1. (i) Assuming (PI), (Strong) is sufficient for either stability of shock waves 
connecting to states along the backward 1-Hugoniot or uniqueness of Riemann solutions involv- 



ing states along the backward 1-Hugoniot; (ii) assuming (P2) ; (Weak) is necessary for stability or 



uniqueness on the backward 1-Hugoniot. (Hi) Assuming (P1)-(P3), (Mediums) is necessary and 



sufficient for stability and (Mediumu) is necessary and sufficient for uniqueness along the backward 
1-Hugoniot through 

(3.1) 



with, moreover, 



Strong 


) => ( 


Mediumu 


) => ( 


Mediums 


) => ( 


Weak 



(iv) Assuming also (P4), there is at most one stability transition as U progresses along the backward 
1-Hugoniot through U+, that is, if the signed Lopatinski determinant turns negative once, then it 
stays negative. 



Proof. Assertion (i) follows from (Strong') and the fact that p > [p] for p_ > 0, assertion (ii) 
from (Weak') in the limit as p- — > 0. In assertion (iii), necessity follows from (Lop)-(Lop a i t 
and (Monotone])^] together with the observation that, by ( |2.16 ), e+ 
e_,p_ — > 0, so that, for e_ > 0, a 2 = — [p]/[r] = p\/e + . 



[t] in the limit as 



Sufficiency follows, similarly, from ( Lop )-( |Lop a it ) and (Monotone), and the observation that their 
righthand sides are monotone increasing with respect p_ and monotone decreasing with respect 
to T-0 so that its minimum is achieved at p- = e_ = and r_ = r max := 2e+/p+ + r+Jj the 



implications (3.1) follow readily from Lax condition |<r| < c. If also (P4) holds, then the righthand 
sides of all conditions are decreasing along the backward 1-Hugoniot curve, yielding (iv). □ 

3.2. The assumptions of Smith. We now recall the assumptions of Smith and show that they 
imply (P1)-(P4). These include, for r > and the structural conditions: 



(Gl) 
(G2) 
(G3) 
(G4) 
(G5) 
(G6) 

Assuming (G6), we may solve p 



e(r, S) > 0. (Positive energy) 
p = — e T (r, S) > 0. (Positive pressure) 
T = e s (r, S) > 0. (Positive temperature) 
p T = — e Tr (r, S) < 0. (Hyperbolicity) 
p TT = — e TTr (r, S) > 0. (Genuine nonlinearity) 
ps = —E t s(t,S) > 0. (Weyl condition) 



-e T (r,S) for S = S(r,p), to obtain e = e(r,p) := e(r, S(r,p)). 



oo, 



Besides the above structural assumptions, Smith imposes the asymptotic conditions: 
(HI) lim e(s,r) = 0, lim e(s,r) = 

s— >— oo s— >+oo 

(H2) lim — e r (s,r) = 0, lim — e T (s,r) 

s— >— oo s— >+oo 

We shall not require Smith's further assumptions that 
(H3) lim T _ 5>0 + e(r,p) = oo for p > 0, 



oo, 



Here, we are using the standard fact that uniqueness of Riemann solutions is equivalent to monotonicity of the 
1-Hugoniot curve (hence, by symmetry, the 3-Hugoniot) in terms of pressure vs. velocity [51 ISml |MP| . 
Note that dependence on (r_,p_) is through factors cr 2 /[p] = — l/[r], ff/[p] = 1/ \J— [p] [r], and l/[p]. 



Here we are using that, by (2.161, t_ is monotone decreasing with respect to both p_ and e_, so that its 



maximum is achieved at the minimum (p_,e_) = (0,0). 
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and 



(H4) 



lim e(s, t) 



0, 



which appear to be used for existence and not uniqueness in [5]. 

We complete our treatment with the following streamlined version of observations of ]W\ E] . 



Proposition 3.2 ([23 EI)- Assumptions pi]) -flG6|) , (pgl|>-pB|| imply fPl 
(i) the backward 1-Hugoniot is C , extending till S — > — oo 
and p, S, and a 2 decrease along the backward Hugoniot. 



-(P4); in particular, 
p — > 0, e — > 0, (ii) t and v increase 







Proof. Recall that the Lax condition is c 2 _ < a 2 < c\ (with the additional condition a < 
distinguishing 1-shocks from 3-shocks). By standard hyperbolic theory |Sm| . (G4)-(G5) imply that 
in the vicinity of (r+, S + ), the backward 1-Hugoniot set of states (t_, SL) satisfying (2.16) and the 
Lax condition consists of a smooth curve on which [t] < 0, [S] > 0. Let us focus on the region 
[r] < 0, therefore. Rewriting the Hugoniot relation as e + — e = (l/2)(p + p + )[r], we see readily 
that, assuming (|Gl|)-(|G2"|), t_ < r max := 2e + /p + + r+. whence r + < r_ < T max . 

We shall show that the backward Hugoniot curve in the vicinity of (r+ , S+ ) extends as a smooth 
curve, monotone in r and S, globally on r G [r + , r max ), 5 € — oo) on which (r, 5) satisfy both 
(2.16) and the Lax condition. Moreover, we shall show that this curve contains all of the points 
satisfying ( 2.16| ) (without the Lax condition) on the region [r] < 0. From these two latter properties, 
we find, reversing the roles of (t-,S~) and (r + ,5+), that any points on [r] > satisfying (2.16) 
consist of reversed Lax shocks, so are not in the backward Hugoniot set, and so we can conclude 
that the curve we have constructed is the entire backward Hugoniot set. 

Applying (|G5|)-p6|), we find that 

[p] = {p( t +,S+) -P{i~+,S-)) + (p(t+,5_) -p(t-,S-)) >Pt(t-,S-)[t] 

whenever [S] > 0, whence — [p]/[t] > p r (r_, S-) whenever [r] < 0, or a 2 > c 2 _, verifying the first 
Lax condition so long as [S] > 0, [t] < 0. 

(i) Differentiating (2.16), we have H$_ = —es_ + \ps- M < by (G3) and (G6), whence, by the 
Implicit Function Theorem, the local backward Hugoniot curve extends as a graph S{t) so long as 
S remains finite. Continuing, we compute that 



+ \(Pr.[r] 



(p+ + P- 



)) = \(p, 



< 



dS- 



< 



so long as [S] > 0, by the (already- verified) first Lax condition, c_ < a A , and thus — —jj^- 

so long as [S] > 0, from which we find that sgn[5'] = — sgn["r] > persists along the entire curve, 
i.e., so long as S remains finite. Rewriting the Hugoniot relation as eo — e = (l/2)(p + Po)[ r ]) we 
see readily that, assuming (G1)-(G2), 

(3.2) r_ < r ma . 



2e /po + to- 



whence r + < r_ < r max , with S — > — oo as r — > r max . Thus, the backward Hugoniot consists of the 
points (r, S(t)) for r G [r + , r max ) and this curve terminates at (t, S) — > (T max , — oo) with e, p — >• 0, 
by (HI ) (|H2h . Moreover, the global inequality H$_ < implies that, as claimed, there are no other 



< 



solutions of ( |2.16 ) for [r] < 0. 

(ii) We have already shown that dS/dr < 0, whence we obtain also dp/dr = ps(dS/dr) + p. 
0, by p T < and p s > 0. This verifies monotonicity of r, p, and S. From [v] = —\/—[p] [t] 
(computed from ( RH[ ) ) , we obtain also the stated monotonicity of v. Next, we compute that 



d 

dS- 



[P] 



[ T ](d/dS-)]pnp](d/dS-)[T] 



has the sign of 



[r](d/dS.)\p]-\p}(d/dS.)[r} = ([p] 



■Pr[T])(dr/dS-) - [r] Ps _ = [r]((c 2 _-a 2 )(dr/dS) 
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PS) 



or, substituting (dr/dS) = -H S _/H T _, [r]( ^_^j ( es +^ s[r]) - ps ^ = -2e s < 0. Thus (since S 

is decreasing) a 2 = — [p]/[t] is monotone decreasing along the backward Hugoniot. Noting that the 
initial value of a 2 at (t + , S+) is c+, we have that c+ > <5 2 , verifying the second Lax condition along 
the backward Hugoniot and completing the proof. 

□ 

Remark 3.3 (Monotonicity of S, a along forward Hugoniot). Symmetric computations show that 
H T+ = \{a 2 - c\)[t] > for [r] < 0, where we are here using the fact already shown that the Lax 
conditions hold whenever [r] < 0, which yields that S parametrizes the forward Hugoniot curve. 
Similarly, (d/dS + )^f\ = —2es + < 0, so that S and a 2 are increasing along the forward Hugoniot 



curve, recovering the results of Weyl [Wj under assumptions (G1)-(G6). 



Remark 3.4 (Monotonicity of t along forward Hugoniot). Monotonicity of r does not always hold 



in the forward direction. Differentiating (2.16) with respect to + variables, we have 



H T+ = e T+ + ^( Pt+ [t] + (p+ +p_)) = \{p T+ [r] - [p]) = -[r](4 - a 2 ) < 

by the second Lax condition, a 2 < c 2 ,, and H$ + = es + + \ps+ [ T ]> hence r is a monotone decreasing 
function of S along the forward Hugoniot curve precisely so long as —H s + /H T+ < 0, or ps < 



-2T/[t\. Recall that this is Erpenbeck's condition, mentioned in Remark 2.14, which is already 



sufficient to imply (Lop 1 ) and stability. Likewise, v and p need not be monotonically related. 



3.3. General pressure laws. We now consider cases of general pressure laws defined in the ab- 
sence of a complete equation of state e. We look at two types of such laws. 

3.3.1. We begin with the question of when a general pressure law 

p = p(r,e) 

determines a full equation of state e = e(r, S), or integral of the fundamental thermodynamic 
relation 

de = TdS — pdr, 

with positive temperature T = e<j, or, equivalently, when we may construct an entropy function 
5 = S(t, e), with S e = e~ l > 0, S(r, •) : M + — > R, satisfying the transport equation 

(3.3) S T -p(T,e)S e = 0. 



As noted in Remark 2.7, for p G C , (3.3) may always be solved locally to particular values 
(ro,eo), by prescribing an arbitrary monotone increasing profile S(tq, •) at To, taking M + — > M, and 
evolving in forward and backward r by the method of characteristics, so long as characteristics 
continue to cover the full halfdine e G M + (they do not intersect, by uniqueness of solutions 
of the characteristic equation de/dr = —p(r,e)), monotonicity being preserved by construction. 
Equivalently, we may prescribe an arbitrary positive temperature profile T(tq, •) at To, such that 
J(l/T)de diverges at both e — > + and e — > +oo. 

The following result, the key observation of this section, states that under the Weyl condition 
we may obtain also global existence in the direction of decreasing r. 

Lemma 3.5. For p G C 1 (M + x M+) satisfying the Weyl condition p e > and p(r, •) = 0, and 
any temperature profile T(tq,-) > such that f(l/T)de diverges at both e — > + and e — > +oo, 



there exists a unique entropy function S = S(t, e) satisfying (3.3) on < r < tq, e G M + , 
with T(r, e) := (S^) -1 > 0. Equivaelently, there is a unique equation of state e(r,S) defined for 
< t < tq and S G M with T = (Se)" 1 = e s and p = —e T . 
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Proof. Under p e > 0, we find that characteristics e(r), de/dr = —p(r,e) diverge in the backward 
(decreasing r) direction, hence continue to cover a half-line [a(r), +00). Moreover, the property 
p(r, •) = ensures that the characteristic emanating from (to,0) remains identically zero, so that, 
by continuous dependence of characteristics, a(r) = and the property lim e _ !> o+ S(t, e) = —00 
is preserved. Meanwhile, the property lim e ^. +00 S(r, e) = +00 is preserved by monotonicity of 
characteristics with respect to e and by the fact that S by construction is constant on characteristics. 

□ 



Corollary 3.6. For a C pressure law p, assumptions (J1)-(J4) imply (P1)-(P4) 



Proof. Applying Lemma 3.5 for any fixed right state r+, we obtain an equation of state e satisfying 
(G1)-(G6) an (H1)-(H2) for < r < 

Tmaxj where T max is defined as in (|3.2l), which generates the 
pressure law p through the relation p := — e T . Noting that only t+ < t < r max is relevant for the 
backward Hugoniot curve through U+, we find by Proposition 3.2 that this implies (P1)-(P4). □ 



Remark 3.7. Under the slightly extended set of assumptions (G1)-(G6), (H1)-(H4), Smith [S] 
obtains not only uniqueness but also existence. It would be interesting to determine what conditions 



beyond (J1)-(J4) would be necessary to obtain existence for a general pressure law p. 
3.3.2. Secondly, we briefly consider the situation of a pressure relation 

p = p(r, T) 

defined in terms of temperature T instead of internal energy e, for example an ideal gas law 

p = p(r, T) := RT/t, R > constant. 
Using a result of [CFJ on the form of possible extensions of the pressure law to a complete equation 



of state, (see (6.77)-(6.80) of [CFJ), Smith [SJ shows that an ideal gas necessarily satisfies (Strong), 
hence the form of the pressure law imposes both uniqueness and shock stability, independent of the 
particular extension to a full equation of state e = e(r, S). 

We note that similar computations may be carried out for an arbitrary relation p = p(r,T); 
namely, conversion to a general pressure law p = p(r, e) amounts to the construction of an energy 
function 

e = e(r,T), e T > 0, 

from which we obtain T = T(r, e) by inversion with respect to e, and thus p(r, e) := p(r, T(t, e)). 
This yields the relations 



(3.4) p T =p T +p T T T , 
Recalling the thermodynamic relation S T 



p e =p T T e , T T = - 
pS e = and S e 



eT 
= f- 



f 

eT 

we obtain T T 



pT e + Tp e = 0, 



hence, substituting from (3.4), 

(3.5) e T = Tp T - p, 

hence e~Tr = &tT = TpTT, and we can solve globally for a positive temperature e and derivative 
for t > tq provided Ptt > 0, and p(r, 0) = 0, which together imply also TpT — p > 0. For further 
discussion, particularly of relation ( |3.5| ), see the discussion of Helmholtz energy in Appendix [Bj 

In the ideal gas case, we have equality, Ptt = Tpr — p = 0, so find that e(r, T) = g(T) depends 
on T alone, where g is any monotone increasing function one-to-one on M + to M + . We thus have 
T T = 0, and, from (3.4), p T = p T , so that (Strong) is equivalent to p T < 0, independent of the 
choice of g. Moreover, we have by direct computation p T = —RT/t 2 < 0, recovering somewhat 
more simply the result of Smith. 

However, note that any admissible temperature function e(r, T) may be modified by the addition 
of any nondecreasing function h(T) with h(0) = to yield another admissible temperature function 
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ei(r, T) := e(r, T) + h(T). Moreover, except in the special situation e T = just considered, or when 
px = 0, such modification will result in a change in T T and p T , so that satisfaction of our conditions 
( Weak ) , ( Strong ) , ( Mediums ) , and ( Mediumu ) may be affected by the choice of e, in contrast with 



the situation of pressure laws of "primary" type p = p(r, e 

Tp 2 

Lemma 3.8. In terms of a pressure law p(r,T), hyperbolicity becomes p T < -g^r, or 



(3.6) p T =p T - p T e T /e T < pp T /e T = pp e 



while (Strong) and (Weak) become, respectively 



(3.7) p T = Pt — Pt^t/^t < (strong) and p T = p T — pt^tI^t < PPt/2&t = PPe/2 (weak). 



Proof. Immediate, from (3.4), □ 

Remark 3.9. These findings are useful for the general van der Waals gas law (p+a/r 2 )(r — b) = RT, 

b 



a,b,R > 0, where p(r, T) = — %. This will be pursued elsewhere. 



3.4. Link to the Riemann Problem. The previous subsections have revealed a certain parallel 
between the transition from stability to instability for individual shocks on the hand and that from 
uniqueness to non-uniqueness for the Riemann problem on the other. Clearly, these two transitions 
are not generally equivalent to each other: there are equations of state for which all shocks have 
positive signed Lopatinski determinant and yet for certain data, the Riemann problem has more 
than one solution. On the other hand, according to Theorem 1.1, the mere existence of shocks with 
negative signed Lopatinski determinant readily implies the existence of Riemann data admitting 
several solutions. This section serves to make this connection transparent. 

As pointed out already by C. Gardner [G] in the context of gas dynamics, vanishing of the 
Lopatinski determinant has a natural interpretation in terms of bifurcation of Riemann solutions. 

Proof of Theorem 1.3. For single-shock data, the Lopatinski determinant coincides with the 
determinant of the local Jacobian of Lax's wave map [Lj. The assertion now follows from the 
following general fact that we briefly prove as we know no reference for it. 

Lemma 3.10. For any C 1 mapping F : MP — > W 1 , local uniqueness prevents det(DF) from under- 
going a local strict change of sign. 

Proof. Local uniqueness at a point p implies the existence of an open ball B around p such that 
F(B) n F(dB) = 0. With d denoting the Brouwer degree, the function B — > R, q h-> d(F, B, F(q)) 
is thus a constant. As d(F, B, F(q)) = sgn detDF(q) whenever detDF(q) ^ 0, B cannot simulta- 
neously contain points q^ with ±detDF(q ziz ) > 0. □ 

Remark 3.11 (Stability as selection principle). The above-described bifurcation analysis bears on 
the question posed in jMPj whether shock stability could serve in the situation of multiple Riemann 
solutions as a selection principle. For, generically, the bifurcation local to the single-shock solution 
should be of fold type, with A taking different signs (corresponding to opposite orientations) on 
either side of the fold; see section 6.2, |Zlj . for an analysis of the general 2x2 case, and |B-G] for 
an extension to the 3x3 gas dynamical case. Assuming that additional small-amplitude waves in 
other family are stable, as is expected to be the case, we find that sign of A, and thus stability of the 
large component shock near (C7_, U + ,a), determines stability of the associated Riemann patterns, 
so that, at the level of Riemann patterns, this is a transcritical bifurcation featuring exchange of 
stability from one branch to another, and indeed serves as a selection principle. 

19 



4. Gas-dynamical examples and counterexamples 



4.1. Global counterexample. Recall the equation of state 

(4.1) E (r,S) = - + C 2 e s / c2 - T / c \ C»l 

T 

asserted to exhibit instability in Theorem A. 

4.1.1. Proof of Theorem 1.2. The function e is evidently convex. Computing, we have 

T = e s = e s /r + e s / c2 - T / c >0, 

p = -e T = e s /r 2 + Ce s / c2 -^ c > 0, 
ps = -e Sr = e s /r 2 + C^eV^/C < 0> 

p T = -e TT = -2e s /t 3 - r 



.S/_:i „S/C 2 -t/C <q 



Pt 



6e 5 /r 4 + C -l e S/^-./C7 >0) 



whence this model satisfies the basic properties (G1)-(G6) assumed by Smith: positivity of e, T, 
p, genuine nonlinearity, and the Weyl condition. 
Continuing, it is readily seen that 



lim e(S, t) 

S->— oo 



0. 



lim e(S, t) 



oo, and lim e(S, r) 



oo, lim e(S, r) = 0, 



validating Smith's asymptotic hypotheses (HI) and (H2). Boundedness of p as r — > + clearly im- 
-oo, or else the first term /r 2 would blow up. Noting that for e s /t 2 < C, e s jr < 



plies that S 1 - 
Cr — )• as t 



H 



we thus obtain e — >• 0, verifying Smith's hypothesis (H3): lim T 
Likewise, (H4) is readily verified: lim T ^ +oc e(s, r) 



0. 



0. Thus, e satisfies all of Smith's hypotheses 



(|G1J-(|G6|) and ([H1JHH4J). 

By Proposition 3.2 and Remark 3.4, we thus obtain monotonicity of S and a along forward and 
backward Hugoniot curves. Evaluating at r + = 1, S+ = 0, and using the assumed asymptotics 
C » 1, we obtain 



-esv 



2e 5 = l + Q(C- 1 ) 2 + Q{C- 1 ) 1 
-e T 3 + OiC- 1 ) C + 0(1) ~3 



> 0, 



giving failure of Smith's condition (Weak). 



4.1.2. Numerical illustration In Figures [T] and [2] below, we illustrate the theory presented in 
Theorem A by numerical computation, displaying the forward and backward Hugoniot curves for an 
unstable 1-shock of (Q, with endstates (r+, S + ) = (1, 0) and (t_,5_) = (12.089548257354499, -6), 
described as plots of r, p, v, e, a against S. Here, the Hugoniot curve is obtained numerically by 
solving (2.16) for r as a function of S using the bisection method. The transition to instability is 
computed along the backward Hugoniot using ( Mediums ) . 



4.2. Local counterexample. Taylor expanding the equation of state (4.1) in the quantitiy 1/C 
about the limiting value l/C = 0asC— > oo, we obtain the simplified, local model 



(4.2) 



e(r, S) = e s /r + S - Ct + t 2 /2, C » 1, 



described in Remark 1.5 It is readily seen that the flow of (2.8) is independent of the value of C, 
with C only serving to shift pressure by positive C but not affecting dynamics. Thus, dropping the 

d energ; 

e(r, S) 



assumptions of positive pressure and energy, we may replace (4.2) by the canonical model 
(4.3) 



e s /r + S + t 2 /2 
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Figure 1. (a) The forward 1-Hugoniot curve through (t_,iSL) « (12.09,-6) for 
global model e(r, S) = e s /t + e 5 ^ _r / c ' of points (r, S 1 ) to which (r_, 5_) connects 
by a Lax 1-shock, displayed as a graph (r,p,v,e,a) over S 1 plotted with respective 
colors (black, green, blue, red, cyan). We zoom in to see (b) the Hugoniot curve, 
(c) T over S, (d) p over S, (e) a over S, and (f) i; over S. Note nonmonotonicity of 
t and v with respect to S. 




Figure 2. The backward 1-Hugoniot curve through (t + ,S+) = (1,0) for global 
model e(r, S) = e s /r + e s ^ c ~ T / C of points (r, £) connecting to (r + , 5+) by a Lax 1- 
shock, displayed as a graph (r,p, v, e, a) over S plotted with respective colors (black, 
green, blue, red, cyan). We zoom in to see (b) the Hugoniot curve, (c) T over S, (d) 
p over S, (e) a over S, and (f) v over S\ The value of (r_, S 1 -) along the backward 
Hugoniot curve at which transition to instability occurs is marked by a black X. 



for which C no longer appears, capturing the large-C behavior of ( 1.3 ). This is perhaps the simplest 
model exhibiting a convex entropy, monotonicity of shock speed and S along the forward Hugoniot, 
but also unstable Lax shock waves. 

Model dOb satisfies dGlb and dG3b-dG6b, but not dG2b, pil-pl, or dPlb-<[P3b. From (pl- 



(G6), we find by the proof of Proposition 3.2 that the backward Hugoniot through any (t+,S+) 
extends for all S G [S+, — oo) as a graph of r as a function of S, on which the Lax conditions hold 
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Figure 3. Extended plots of the (a) forward and (b) backward Hugoniot curves 
shown in Figures 111 and v2j for global model e(r, S) = e s jr + e s l c ~ T / C . 



and p, a 2 are monotone decreasing; likewise, the forward Hugoniot is parametrized by increasing S, 

with a monotone decreasing (a 2 increasing). However, p and e do not approach zero as S — > — oo, 

instead becoming negative. 

Indeed, the Hugoniot curve for (4.2)-(4.3) may be solved explicitly as a cubic in r. For example, 

s 

setting (r+, S+) = (1, 0) and noting that p(r, S) = —e T = ^-r, we find that p + = and e + = 3/2. 

Substituting into H(t, s) = e — e+ + (l/2)(p + p + )(r — r+) = yields 



H(t, S) 



or, multiplying by r 2 , 
(4.4) 



+ S' + l/2r 2 -3/2 + l/2 



r 3 + 3e s r + (25-3)r 2 - e 5 



r ( T _ 1) = 



0. 



Thus an explicit solution t(S) is available through the cubic formula. Moreover, along the backward 
Hugoniot, we find as 5 — > — oo that (4.4) reduces asymptotically to r 3 + 2SV 2 = 0, so that, 
asymptotically, r(S) ~ —25 — > +oo, and e ~ 2S 2 — > +oo, p ~ — r 



-oo. 



Proposition 4.1. For model (4.3) for, equivalently, (4.2)), the backward Hugoniot through any 
(t + , S+) extends for all S G [S 1 .)-, — oo) as a graph of t as a function of S, on which the Lax condi- 
tions hold, r is monotone increasing, and p, a 2 are monotone decreasing; likewise, a is monotone 
decreasing and S monotone increasing along the forward Hugoniot. Yet, for (r+,5+) = (1,0), 
there exist (t_,5L) along the backward Hugoniot such that the Lax 1-shock connecting (r,S)± is 
Lopatinski (i.e., Hadamard) unstable. 

Proof. By direct computation (see further calculations just below), we obtain 

-e ST _ 2es\ 1 4 1 



as S 



-oo, r 



+oo, and p- 



[p] J (t,S) = (t+,S+) 



> 0, 



-P- 



(Weak'). 



-oo, yielding the result by Proposition 2.10 and failure of 

□ 



4.2.1. Estimating the stability transition. From (4.3), we obtain 

p = e s /T 2 -T, e s = e s /T + l, e r 
hence, evaluating at (5+,r+) = (0, 1), 

P+ = 0, (es)+ = 2, (e rr )+ = 3, (e r5 ) 
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2e s /r 3 + 1, e rS = -e s /r 2 , 



so that 6 := 



_ (2)(3) 

— k oti^ a — 1 _i_ 



f Z = 6,and0:=l + V^ = l + V ^ IJ . 



Thus, change in stability correponds to change in sign of 



i:=p + -p_-^ = ~+T-6-2 1 /- - 3 r ) (r- l.) 




Noting that p ~ — r, for 5- large and negative, so that [p]/[r] ~ - , and recalling that c: 



2 



—p Tt + = 3 we may well approximate ~ 1 + y 6 ( r _i ) ~ 1-6, so that the stability transition occurs 

approximately at r_ = 9.6 and, estimating 5 « 4(3 — r) for e s << 1 in (4.4), 5_ rj —3.3. In 
fact, numerically we determine that 5 = occurs for 5_ G (51,52) = (—3.3348293,-3.334829), 
r_ € (t"2,ti) = (9.6589763,9.6589769). See Figure [4] below for a plot of the Hugoniot curve 
computed by the Bisection Method, with the stability transition indicated on the curve. 

Remark 4.2. From c\ = 3, a 2 = — M ~ - T Z T+ ~ 1 + 7-, and c 2 _ = 1 + 0(e s ~), we verify directly 
the Lax conditions (?_ < a 2 < A. for 5- » 1, with the shock nearly characteristic at U-. 
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(a) 



(b) 



Figure 4. We mark transition to instability with a thick X. (a) Plot of the Hugoniot 
curve for the local model, e(5,r) = e s /t + S + t 2 /2, (r ,5 ) = (1,0). (b) Plot of the 
the Hugoniot curve for the local model with a thick red line and that of the global 
model for C = 40, 100, 250. Note that as C — > 00, the Hugoniot curve of the local 
model matches well that of the global model. 



4.3. Stable example. For comparison, consider finallly the closely related model 



(4.5) 



e(r,5) = e s /T-Cr + T 2 /2, C » 1. 



Here, likewise the flow of (2.8) is independent of the value of C, so that it is equivalent to consider 
the simpler version 

(4.6) e(r, s) = e s /r + t 2 /2. 



Proposition 4.3. For model (4.6) (or, equivalently, ( |4.5| ) ), the Hugoniot through any (t+,5+) 
extends for all S G (—00, —00) as a monotone graph of r as a function of S, on which the Lax 
conditions hold precisely for [r] < 0. Moreover, all Lax shocks have positive signed Lopatinski 
determinant. 



Proof. Again, from (G3)-(G6), we find by the proof of Proposition 3.2 that the backward Hugoniot 
through any (r+, 5+) extends for all 5 G [5 + , —00) as a graph of r as a function of 5, on which the 
Lax conditions hold and p, a 2 are monotone decreasing; by similar reasoning, the forward Hugoniot 
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is parametrized by increasing S, with a monotone decreasing (a 2 increasing). Indeed, the Hugoniot 
curve for (4.6) can be determined explicitly by solving 

= H(r, S) = e s /r + r 2 /2 - e+ - \ (e s /r 2 - t + P+ )(t + -t) 



e s /r 



\{e s /t 2 + P+ ){t + -t) + \tt + 



or, using the identities p + + t 4 



£± and e+ + p -±f± = |^±, e 3 ^ 

t+ ' 1 z r_)_ 7 



3t_|_— t 
3t— t_|_ 



for 



S-S+ = 

6(x-l) 2 



log 



3r 4 



T\ T 



3r — r 4 



^2 ■ 



Computing (d/dT)(S — S + ) = — x ^°x)(3x-i)t + < f° r x := T / r + an d 1/3 < x < 3, we find that 
S is a monotone function of r along the entire (forward and backward) Hugoniot curve. But, 
monotonicity along the forward Hugoniot implies stability by Remarks |2.14| and 3.4 □ 




Figure 5. (a) Plot of the Hugoniot curve for the local model e(S,r) = e /r — 
Ct + t 2 /2 with C = 100, (to, Sb) = (1,0). (b) Plot of 5 against S where 5 < 
corresponds to stability. 



PART II. VISCOUS STABILITY 

Stepping outside the inviscid setting, we can extract further interesting information from the 
sign of the signed Lopatinski determinant. This fact and the Evans function, to which it refers, 
are discussed in Section 5. Section 6 presents the findings that were summarized as Numerical 
Observation A in Sec. 1. The phenomenologically most interesting piece of Part II might be Section 
7, where the findings of Numerical Observation B are detailed, on the rich stability behavior of 
shock waves in an artifically designed system. 

5. The Evans function 
We now turn to the study of viscous stability, or stability of continuous traveling shock fronts 

(5.1) w(x, t) = w(x — at), lim w(z) = U± 

z— >±oo 

of a hyperbolic-parabolic system of conservation laws 

(5.2) f(w) t + f(w) x = (B(w)w x ) x , 
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corresponding to an inviscid shock (U-,U+,a) of (2.1), that is, a solution of the traveling-wave 
ODE 



(5.3) B(w)w' = f(w) - af°(w) - (/([/_) - vf(U-)). 

Recall |GZ1 \ZH\ IMaZH lMaZ2l |Z1 | IBLeZl |Z5] that, under rather general circumstances, shared by 
most equations of continuum mechanics- in particular, the equations of gas dynamics, MHD, and 
viscoelasticity- Lax- type traveling- wave profiles w, when they exist, (i) are generically smooth and 
converging at exponential rate as z — > ±00 to their endstates U±, and (ii) are nonlinearly orbitally 



stable as solutions of (5.2) if an only if they satisfy an Evans function condition: 



(5.4) has no zeros on KA > except for a multiplicity one zero at A = 0, 

where the Evans function D is a Wronskian at x = of appropriately chosen bases of solutions of the 
linearized eigenvalue equations about w that decay at plus and minus spatial infinity, respectively, 
that is analytic in A on a neighborhood of 3?A > 0, with complex symmetry D(\) = D(A), whose 
zeros agree on {3ftA > 0} \ {0} in location and multiplicity with the eigenvalues of the linearized 
operator about the wave. 

5.1. Viscous vs. inviscid stability. Like inviscid stability, viscous stability holds always for 
Lax type shocks in the small-amplitude limit \ZH\ IHuZ3|, IPZl IFS1} IFS2] . Moreover, as found in 
|GZl IZS| , we have the fundamental relation 

(5.5) D{\) = Xv5 + 0{\ 2 ) 

for A sufficiently small, where v is a Wronskian associated with the linearized traveling-wave ODE, 
with nonvanishing of v corresponding to transversality of the connecting viscous profile, and 5 



is the Lopatinski determinant (2.4) associated with the inviscid stability problem. Without loss 
of generality, normalize v and 5 so as to be real. Then, moving along the Hugoniot curve in the 
direction of increasing amplitude, we find that, so long as there persists a transversal traveling- wave 
connection, so that v cannot change sign, that real zeros of the Evans function can cross through 
the origin into the unstable half-plane if and only if the Lopatinski determinant 5 changes sign, 
that is, there is a corresponding stability transition at the inviscid level |GZ1 \ZH\ IZSl \Z1\ IZ3j . 

More precisely, we see that, assuming persistence of transversal connections, inviscid stability 
transitions correspond at viscous level to the passage of an odd number of unstable roots through 
the origin, so can never occur without a corresponding viscous stability transition. On the other 
hand, there can occur viscous stability transitions corresponding to passage of an even number 
of real roots, or of one or more complex conjugate pairs of unstable roots through the imaginary 
axis, which are not detected at the inviscid level. Thus, viscous stability is a logically stronger 
condition than inviscid stability [ZSJ. Our goal in this part of the paper is to investigate whether 



this logical implication can be strict for systems (5.2) possessing an associated convex entropy, 
and in particular whether complex conjugate pairs of roots can cross, signaling Hopf bifurcation to 
time-periodic "galloping" behavior [TZTl [TZ2l [TZ3l iBeSZ] . 

5.2. Numerical stability analysis. The Evans function, or, for that matter, the underlying 



profile (5.1), is only rarely computable analytically. However, as described, e.g., in IBrZ|,IHuZH 
Z4j . these may be approximated numerically in a well-conditioned and efficient manner. Here, we 
carry out numerical experiments using MATLAB and the MATLAB-based Evans function package 
STABLAB [BHZ2J developed for this purpose by Barker, Humpherys, Zumbrun, and collaborators. 
More precisely, we study the integrated Evans function -D(A), a Wronskian associated with the 
integrated eigenvalue equations, whose zeros on {3ftA > 0} \ {0} likewise agree in location and 
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multiplicity with the eigenvalues of the linearized operator about the wave, but which does not 
have a zero at the origin, obeying a small-frequency expansion 

(5.6) D(X) = u5 + 0(X) 



in place of (5.5) 



Descriptions of the underlying algorithms and other information useful in reproducing these 
computations may be found in Appendix [Cj 

6. Gas dynamical examples 



The compressible Navier-Stokes equations, augmenting (2.8) with the transport effects of vis- 



cosity and heat conduction, appear in Lagrangian coordinates as 

n-v x = 0, 

(6.1) v t+p x = (—) x , 

rHVV x 



(e + v 2 /2) t + (vp) x 



), + (^) s 



T T 

where r denotes specific volume, v velocity, e specific internal energy, p pressure, and T temperature, 
and n and k are coefficients of viscosity and heat conductivity, here taken constant. 

To close these equations requires both a pressure law p = p(r, e) and a temperature law T = 



T(t, e), in contrast to the case of the Euler equations (2.8), that is, a complete equation of state. 
Indeed, the most convenient formulation for our purposes here is in terms of the variable w := 
(r, v, T) replacing e with T. We will use equations of state given in the form 

(6.2) p = p(r,T), e = e(r,T) 

arising through the Helmholtz energy formulation described in Appendix [B| for which the equations 
take the convenient block structure described in Appendix |C.1[ Note that besides hyperbolicity, 
p T < 0, the Navier-Stokes equations (6.1) require also parabolicity, 

ess = (dT/de)(de/dS) > 0, 



in order to be well-posed, which is the condition needed to obtain laws (6.2) from a given energy 
function e = e(r, S) of the form considered in previous sections. 



Substituting (5.1) into (5.2) and solving the first equation of (5.3) for 
(6.3) 



u>i 



a 



as described in Appendix |C.l[ we obtain for a profile corresponding to an inviscid shock (Z7_, U+,a) 
of (2.8) the traveling- wave equations 



(6.4) 



-<X[V 



-a(e + v 2 /2 



) + p- p- 
vt/2) +pv 



P-V- 



Proposition 6.1 (|Gi|). For a C 2 energy function e(r,S) satisfying (G3)-(G6) ; every Lax type 
shock of (2.8) has a unique continuous shock profile (5.1) of (6.1)-(6.2), which is a transversal 



connection of (6.4) converging exponentially to endstates U± as x — )• ±oo in up to two derivatives. 



Proof. More generally, Gilbarg shows that, assuming the Weyl condition ps > 0, a Lax-type inviscid 
1-shock for which the pressure at the righthand state U+ is minimal among all connections to U- 



with the same shock speed a possesses a unique viscous profile (5.1). Since (G3)-(G6) imply 



uniqueness of inviscid Lax shocks for a given speed, we thus obtain existence of a unique profile; 
transversality then follows by the observation that this connection must be of saddle-node type, 
so transversal whenever it exists and exponentially decaying by standard stable/unstable manifold 
theory. See also [W| and jMP| Appendix C]. □ 
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Proposition 6.2 (Kaw,KSh). For a convex C 2 energy function e(r, S) satisfying (G3)-(G6), (6.1) 
possesses the convex viscosity-compatible entropy S(t, e), globally onr,e> 0. 

Proof See, e.g., [Kawl iKShl IMP! IZl] [Z2] . 



□ 



6.1. Numerical details. Despite the familiarity of the gas dynamical equations (6.1), and the 
straightforward analytical treatment under the assumptions we impose, we find the numerical 
Evans function analysis of our gas dynamical examples to be quite challenging, due to the presence 
of multiple scales. Indeed, this turns out to be the most computationally intensive problem that 
we have so far considered, even including the famously intensive problem of detonation stability 
considered in |Er31 ILSl BHLyZl, IBZ2j . We thus find it necessary to include some extra details on 
how the computations are carried out, despite that this is not our main emphasis here. 

6.1.1. Computation of profiles. To approximate numerically the profiles guaranteed by Proposition 



6.1 for a given equation of state e(Y, S), we see two basic possible approaches. The more direct, but 
numerically somewhat cumbersome approach is to evaluate the functions p(r, T) = p(r, S(t, T) and 
e(r, T) = e(r, S(t,T) by numerically inverting the relation T = e~s(r,S) to solve for S = S(t,T), 
using ess > 0. As computation of the profile is a one-time cost, the associated cost in function 
calls is perhaps not too serious, but this still seems preferable to avoid. 

The approach we follow here is, rather, to compute the profile in (r, S) coordinates using the 

1 ' ' 1 deriving from T' = es T T~' + essS' and the 



esr e S s 



change of variables formula 
relation r' = — a v' obtained from (6.3), yielding, finally, 

-a 

esr e S s 



(6.5) 



//;> 



-a(v 



a(e + v 2 /2-e- 



+ p - p- 
W-/2) + pv — 



From (f, S)(x), we can recover all other variables through (6.3) and the equation of state. 

A final detail is that in the large-amplitude limit S — > — oo, the traveling-wave ODE (6.5 ) features 
a wide separation of scales. Specifically, the linearized equations about endstates U±, 

l 



(6.6) 



e"ss 




1 

ess 



± 



a 2 + c 2 



v(2a 2 



PS 

-oT + psv 



for global model (4.1) in the regime 1 << \S-\ << C considered in the crucial stability transition 

-oo) to (computing |<r| ~ 1, c = — e TT ~ 1, e~s ~ 1 5 e~Sr ~ C , 



regime, is asymptotic at U. 
ess ~ C" 2 ) 



(x 



v-C 




the regime r_ — > +oo, 5_ 
growth rates ~ 1 j [i and r_ e r ~ 



1/fi and C 2 j k differing by two or more 
• is inherited by the 
oo roughly corresponding to C 



which exhibits growth modes with exponential rates 

orders of magnitude. Similar but more extreme behavior is inherited by the local model (4.3) in 

i(3-r)- 
-I 2 /k. 



oo for (4.1 



with 



We find that this gives trouble for the standard method (described further in Appendix C.2) 
of solving for profiles with a boundary-value solver with projective boundary conditions; indeed, 
we are unable to compute by this method beyond approximately the range \C\ < 75 for (4.1) or 
r_ < 19 for 

Above these values, we rely on a more primitive, shooting method, starting at a point distance 
10~ 8 from the saddle U + in the direction of the direction of the stable subspace of (6.6), and 



This is with initial guess given by a simple tanh interpolation. We can go slightly farther by continuing along 
the Hugoniot curve, using previous profile as initial guess; however, the computations are still extremely slow. 
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integrating backward in x the traveling-wave ODE (6.5) using an appropriate stiff numerical ODE 
solver until the solution comes within 10 -8 of the repellor U-. This simple solution gives excellent 
results for essentially arbitrary values of S- . 

6.1.2. Evans computations. In computing the Evans function, we face the same numerical stiffness 
issues that we faced for the profile equation (recall that at freqency A = 0, the eigenvalue equation 
reduces to the linearized traveling- wave ODE). Though it is standard practice for us to compute 
the Evans function by shooting, we find it necessary to carry out this shooting algorithm using a 



stiff ODE solver rather than the standard RK45. For further discussion, see Appendix C.2 



A related issue is that for C > 10 or so in (4.1 ), or for r_ > 6 in (4.3), the standard Evans function 
experiences exponential growth/decay arising from spatial variation in the coefficient matrix of the 
Evans (i.e., first-order eigenvalue) system so great as to go out of numerical range- in this case, 
returning a value that is to numerical approximation identically zero. Similar phenomena have 
been seen to arise in other numerically difficult Evans computations |BZll IBZ2| BHLyZi] , We 



follow the simple solution here (mentioned also in the cited references) of turning off growth; that 
is, computing using the polar coordinate method of |HuZl] . we suppress variation in the radius 
and simply evolve orthonormal bases of the subspaces of decaying solutions at plus and minus 
spatial infinity (the "polar no radial" setting in the STABLAB Evans package). This amounts to 
the continuous orthogonalization method of Drury |Dr| , which preserves winding numbers /zeroes 
of the Evans function, but loses the desirable property of analyticity [Z4]. (It is in fact C°° in A.) 

The same numerical stiffness issues prevent degrade our usual estimates excluding unstable eigen- 
values outside a semicircle of radius R. In this section, therefore, we often compute for a fixed radius 
R without guarantees that there are no unstable eigenvalues of larger modulus. 

Remark 6.3. The presence of the very slowly decaying term e s / c2 - T / c i n the Evans function co- 



efficient arixing in our global counterexample (4.1) degrades the standard convergence estimates 



described in, e.g., |Z4} IZ5| . wich depends on uniform exponential convergence of the coefficient 
matrix of ( C.7| ) to its limits as x — > ±oo. However, uniform- in-C estimates are easily recovered 



similarly as (but much more simply than) in [HLZJ , by the incorporation of an additional "tracking" 
argument in the slowly- varying far- field regime. A similar argument, though we do not carry it out 
here, may be used to rigorously verify convergence as C — > oo of the Evans function for the global 



counterexample (suitably normalized) to the Evans function for the local example of Section 3.1 



6.2. Global counterexample. Returning to the equation of state 
(6.7) e(r,S) = - + C 2 e s / c2 - T / c , C»l 

T 

shown to exhibit inviscid instability in Theorem A, we report the results of our numerical Evans 
function investigations. In our study we examined the parameter set S- = {— 1, — 1.1, — 50}, 
(r + ,5 + ) = (1,0), for C = 10. We solved the profile using Matlab's boundary value solver bvp5c. 
In addition, we computed the Evans function on a smaller radius for (t + ,S+) = (1,0) and several 
values of S- including S- = — 10 5 for C = 100 with R = 4, solving for the profile using Matlab's 
implicit ODE solver, odel5s. 

The results of these investigations were, in all cases, that, as shock amplitude was increased (i.e., 
S— decreased), there was, similarly as in the inviscid stability analysis, a single stability transition 
consisting of a real root of multiplicity one passing through the origin into the unstable complex 
halfplane KA > 0, where it remained thereafter as the only unstable eigenvalue. In particular, there 
were seen no complex unstable eigenvalues at any point in the process. That is, in this case, viscous 
and inviscid stability predictions appear to agree. 

For C = 10, the computed inviscid stability transition occurs approximately at S- = —23.2, 
whereas the computed viscous stability transition occurs for S- 6 (—22.1, —22.2); recall that, by 
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the analytical theory, these two values should agree, as parameters for which a real root passes 
through the origin. Except for this difference in location of the transition point, the results of 
inviscid and viscous numerical stability analyses agreed. The relatively large size of this discrepancy 
we ascribe to the stiffness of the underlying eigenvalue equation both for viscous and inviscid cases; 
compare for example to the excellent agreement found in the investigations of the nonstiff system 
considered in Section 

The study was further challenged by the slow rate at which the Evans func tion converges to 
its high frequency behavior D(\) ~ Cie C2V ^, as described in (C.8), Appendix C.2 In tables [l] 



and [2| we provide data about curve fitting the Evans function with C\e c ' 2 ^ as A — > oo. As seen 
in the tables, the needed radius for convergence for large C exceeds the range where numerical 
computations can be carried out accurately. Consequently, we arbitrarily chose R = 250 for the 
Evans function study. Figure [6] demonstrates typical output from our study. 



6.3. Local counterexample. Next, we turn to the local model (4.3) obtained by expansion 
of the global model in the C — > oo limit, describing the results of a corresponding numeri- 
cal Evans function investigation in that case. In our study we examined the parameter set 
S- = { — 1, —1.1, — 8.1}\{— 3.3}. Numerically, we found that both the inviscid and viscous stabil- 
ity transitions occurred around 5_ = —10/3. For all parameters covered in our study, the viscous 
Evans function had one root or no root according as the shock was inviscid unstable or stable. 
As with the global model, we could not carry out a numerical high frequency study, and so we 
arbitrarily chose R = 250 for the contour radius in our Evans function computations. Figure [7] 
demonstrates typical output from our study. In conclusion, the results of our numerical study 
indicate in this case too that viscous and inviscid stability agree. 



6.4. Stable example. Finally, we briefly discuss the case of example (4.6) seen to be inviscidly 
stable. In our study we examined the parameter set S- £ {—0.1, —0.2, —11.5}. For a few values 
the boundary value solver could not solve the profile equation, returning a "singular Jacobian" 



R 


error 


Ci 


c 2 


2 


0.8726 


0.03233 


2.427 


4 


1.693 


0.002496 


2.996 


8 


2.103 


3.386e-05 


3.639 


16 


6.33 


3.659e-08 


4.281 


32 


12.47 


1.476e-12 


4.816 


64 


17.38 


1.784e-18 


5.108 


128 


4.211 


3.922e-26 


5.171 


256 


0.9991 


7.957e-36 


5.051 


512 


0.9949 


3.594e-48 


4.828 


1024 


0.9994 


8.181e-65 


4.611 


2048 


0.9945 


2.664e-88 


4.456 


4096 


0.9805 


5.719e-122 


4.362 


8192 


0.948 


3.166e-170 


4.312 


1.638e+04 




NaN 


NaN 



Table 1. For C = 10, S_ = —15 in the global model we compare how well C\e 2%/ 
approximates -D(A) for |A| = R. Here R is the radius of the semicircle on which we 
compute the Evans function, C\ and Ci are the curve fit returned, and the conver- 
gence error is the maximum relative error between D(\) and C\e C2 ^ evaluated at 
X = R, Ri. 
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Figure 6. Viscous study results for e(S,r) = e s jr + C 2 e c ~ 1( - C ls ~ T ^ with S- = 
—30, C = 10, fj, = k = 1, (to, So) = (1)0)- ( a ) Viscous profile, (b) Integrated 
Evans function computed with the adjoint polar coordinate method, evaluated on a 
semi-circle with radius R = 1. Winding number is 1. (c) As in (b) but with R = 250. 
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Figure 7. Viscous study results for e(S,r) = e 5 /t + S+t 2 /2, 5- = -5, /u = « = 1, 
(to, So) = (1)0). (a) Viscous profile, (b) Integrated Evans function computed with 
the adjoint polar coordinate method, evaluated on a semi-circle with radius R = 1. 
Winding number is 1. (c) As in (b) but with R = 250. 



error. For these, we solved the profile using MATLAB's implicit ODE solver, odel5s. To ensure 
that we computed the Evans function on a contour enclosing any possible unstable eigenvalues, we 
performed a high-frequency study curve- fitting with C\e C2 ^, requiring that relative error be no 



R 


error 


Cl 


c 2 


2 


0.3373 


0.711 


0.2412 


4 


0.3457 


0.9109 


0.04668 


8 


0.287 


1.268 


-0.08408 


16 


0.2102 


1.876 


-0.1572 


32 


0.1479 


3.003 


-0.1944 


64 


0.1055 


5.439 


-0.2117 


128 


0.0779 


12.07 


-0.2202 


256 


0.05837 


36.77 


-0.2253 



Table 2. For C = 2 in the global model we compare how well C\e 2 * approxi- 
mates D(\) for |A| = R. Here R is the radius of the semicircle on which we compute 
the Evans function, C\ and C 2 are the curve fit returned, and the convergence error 
is the maximum relative error between D(X) and C\e C2 ^ evaluated at A = R, Ri. 
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greater than 0.2 between the Evans function and the approximating function. Table [3] demonstrates 
high-frequency convergence for a typical parameter. In addition we used an adaptive A mesh to 
ensure that the change in argument was no greater than 0.2 between successive A points. We 
found for all parameters studied that the Evans function has no root in the right half complex 
plane, consistent with stability. The maximum relative change in A for each contour computed was 
0.1005. The average value of R was 4 and the average number of points computed was 40 (80 points 
after using conjugate symmetry of the Evans function to reflect the values computed on the upper 
half-plane 9 A > to the lower half-plane). The average time to solve a profile was 12.8 seconds and 
the average time to solve the Evans function was 6.76 seconds using 8 cores on a Mac Pro. Evans 
function output for a typical point in parameter space is given in Figure |8| In all cases, we found 
a winding number of zero, indicating viscous stability in agreement with our inviscid findings. 

Evans Function 




x " Re(X) Re(A) 

Figure 8. Viscous study results for e(S,r) = e s ' jr + r 2 /2 with (t ,5o) = (1,0), 
and S- = —5. (a) Viscous profile. (b)Integrated Evans function computed with the 
adjoint polar coordinate method, evaluated on a semi-circle with radius R = 250. 
Winding number is 1. (c) The Evans function evaluated at 8 points along a semi- 
circle of radius R = 256 is plotted with open circles. The best curve fit of the Evans 
function with Cie C2 ^ is plotted with closed dots. The relative error between the 
two is 0.069. 
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error 


Ci 


c 2 


2 


0.04501 


0.6019 


0.359 


4 


0.08243 


0.4554 


0.3933 


8 


0.1244 


0.2953 


0.4312 


16 


0.147 


0.1556 


0.4651 


32 


0.1416 


0.06249 


0.4901 


64 


0.1179 


0.01757 


0.5052 


128 


0.09174 


0.002994 


0.5136 


256 


0.06888 


0.0002471 


0.5191 


512 


0.05042 


7.366e-06 


0.5223 



Table 3. For S- = — 5 in the stable model we compare how well C\e 2 * approxi- 
mates -D(A) for |A| = R. Here R is the radius of the semicircle on which we compute 
the Evans function, C\ and C 2 are the curve fit returned, and the convergence error 
is the maximum relative error between -D(A) and C\e C2 ^ evaluated at A = R, Ri. 
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7. Designer systems 

We conclude by investigating the range of possible behaviors for general systems possessing a 
convex entropy. Let A(x) be a given real symmetric n x n matrix-valued function, and B(x) a given 
real positive definite n x n matrix-valued function, in the sense that m :=\{B + B T ) > 0, each 
converging exponentially as x — > ±00 to limits A± and B±, and consider the eigenvalue problem 

(7.1) Xu = d x (B(x)d x - A{x))u, u e C n , 

resembling those found in the stability analysis of viscous shock waves. This might be considered 
as a model for the type of eigenvalue equation that would arise from a symmetrizable system, in 
particular a system possessing a convex entropy. 

Our first observation is that, by an ingenious construction of Bianchini [BJ, any such system 
can in fact be realized as the eigenvalue problem for a viscous shock wave solution of a system 
of conservation laws with a viscosity-compatible strictly convex entropy. Specifically, consider the 
system 

u t + (A(v)u) x = (B(v)u x ) x , 

(? - 2) A Mf ^ 1 

v t + {-u • A (v)u + -v ) x = v xx , 

where v € R 1 . 

This contains a trivial standing viscous shock solution (u, v) = (0, v(x)), where v{x) := — tanh(a;/2) 
is a shock of Burgers equation vt + (v 2 /2) x = v xx . Linearizing (7.2) about (u,v), we obtain 

™ C0,+(T ")(:)^((T ;)(: 

which decouples into a (stable) linearized Burgers equation in v and an equation 

(7.4) u t + (A(v(x))u) x = (B{v(x))u x ) x , 

that has the associated eigenvalue equation (7.1) provided that we choose A(y) := A{2 arctan(— v)) 
and B(v) := B(2 arctan(— v)). Moreover, it is easily verified that (7.2) possesses a viscosity- 
compatible stricly convex entropy rj(u,v) := \\u\ 2 + \\v\ 2 , with associated entropy flux q(u,v) := 

\{u, A((v)u) + 5M 3 + \v(u,A'{v)u. That is, drjdF = dq and d 2 r] ^\ = ^\ is positive 

/ A(v)u \ 

definite, where F(u,v) = ( 1 j^/^ ) + 1 2 J " means that we may search for stability phe- 
nomena of symmetrizable systems possessing a viscosity-compatible strictly convex entropy, by 
considering the linear problem (7.1), which is free to our specifications. 

7.1. The rotating model. A natural choice in seeking examples of instability is the rotating model 
consisting of (7.2) with 

(7.5) A(v) = R e{v) A m R_ e{v) , B(v) = Id, 

/cos (9 -sin(9\ , (\ 0\ 
where K$ := \ . n n , A m = _ 1 , and 

\sm9 cos8 J ' m \0 _ V 

(7.6) 6(v) = Mttv 

M an arbitrary real number, considering the family of stationary shocks Vj(x) := — 7tanh(7a;/2). 
For, it is readily seen that these are Lax 2-shocks, with a Lopatinski determinant 

'cos — M7T7 — sin(M-7T7) \ 
5 7jM = det | sin -Mttj cos(Mirj) = -27(cos 2 (M7T7) - sin 2 (M7T7)) 
-2 7 / 
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that changes sign from negative to positive and back infinitely often as 7 increases from the small- 
amplitude limit 7 = 0, with the first stability transition occuring at Mir^* = vr/4, or 7* = 1/4M. 
This in passing gives another example of a 3 x 3 system with convex entropy exhibiting inviscidly 
unstable shock waves, 



7.1.1. Evans system. For profile u 
decoupled, u equation (7.4) is 



= 0, v = — 7tanh(7x/2), the integrated Evans system for the 



(7.7) 



W = A(x, \)W, A 







I 

B- l A 



(v(x)). 



Remark 7.1. Setting W = TY, where T := 
(7.8) M = T- 1 AT-T- 1 T' 



Re 

R e 

/ 

A Afj 



we may convert (7.7) to Y' = MY , where 



Mirv'(x) 



J 
J 



and J 



-1 

1 



Involving only exponential functions in x entering through the scalar multiplier 



v' , this seems possibly amenable to exact solution, an interesting direction for further study. 

7.1.2. Transversality of profiles. Again appealing to the triangular form of the equations, we find 
that transversality of profiles is equivalent to nonexistence of asymptotically decaying solutions of 
the u-component u' = A(v(x))u of the linearized standing-wave equation. More, the transversality 
coefficient v, like the Lopatinski determinant, factors into —27 times a transversality coefficient for 
this decoupled u-component. 

Lemma 7.2. For \M^\ bounded, profiles are transversal for 7 sufficiently small. More, when 
appropriately normalized, (i) v/{— 27) — > 1 as 7 — > + and (ii) v/(—2j) — > 5/(— 27) as 7 — > +00, 
while (Hi) v/(—2j) changes sign at least [4(M7 — 1)] times as 7 increases from + to +00, where 
[■] denotes the greatest integer function. 



-1 



Proof, (i) (Basic tracking argument; see [Z5].) Setting u = Rez, similarly as in Remark 7.1 we 
obtain 

z' = A m z + Mv'(x)J, with J := ^ 

Observing that |M-u'(x)| < M7 2 < C7 — > for 7 — > and \M~/\ < C, we obtain the result by the 
Tracking Lemma of [ZHJ, a standard estimate for slowly- varying-coefficient systems, after factoring 
out the (real) exponential growth rates in growing and decaying modes. 

(ii) Noting, for (^-Q that Evans system {77} satisfies \A(x)-A±\ < CMje~^ x \ < C 2 e~^ x \ 
for x ^ 0, we obtain by the Convergence Lemma of |PZj that the Wronskian v converges up to an 
error of order H^e -7 ^' ||^i = 0(C/^) — > to the determinant of a (smoothly chosen real) stable 
eigenvector of ^4+ and a (smoothly chosen real) unstable eigenvector of A_, i.e., 5/{— 27). 

(iii) Reviewing the previous arguments more closely, we find that for 7 sufficiently small the 
decaying mode as x — > —00 rotates angle 2Mjtt, while for 7 sufficiently large it rotates at most ir, 
with the directions of stable and unstable subspaces as x — > ±00 held fixed for all 7. Thus, as 7 
goes from to +00, the decaying mode as x — > 00 crosses the decaying mode as x — > +00 at least 
[4(M7 — 1)] times, each time signaling an associated change in sign. □ 



From Lemma 7.2, we find that, fixing K sufficiently large, taking 7* sufficiently small and 
M := K/j*, and varying < 7 < 7*, the profiles u 7 remain transversal, while the Lopatinski 
determinant changes sign ~ 2K times, signaling (recall ( |5.6| )) 2K passages of an eigenvalue through 
the imaginary axis, hence appearance of up to 2K unstable roots. (Our numerical experiments 
indicate that these roots all cross in one direction, similarly as in standard Sturm-Liouville theory.) 
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On the other hand, fixing M7 = constant, and increasing 7 from + , we find that 5 stays 
constant, while v changes sign at least ~ 4(M7 — 1) times. More, varying M7 within a bounded 
set for 7 sufficiently large, we find that changes of sign in 5 and v may be made to occur arbitrarily 
close together, increasing the chance of collision of roots and subsequent splitting into a complex 
conjugate pair. (Recall that changes of sign of 5 and v signal crossings through the origin of roots of 
the Evans function D.) This simple rule of thumb guides our strategy in searching for complex roots 
and Hopf bifurcation, and indeed gives a reasonable a posteriori fit to the data that we see. Why 
some colliding roots split into the complex plane while others appear to stay real is an interesting 
question to which we do not have an answer. 

7.1.3. High-frequency bounds. We start by bounding the modulus of unstable eigenvalues. 



Lemma 7.3. For system (7.3) with the choices (7.5)-(7.6), there are no unstable eigenvalues with 
modulus |A| > 4. 



Proof. Let A be as in 7.5 and consider the integrated eigenvalue problem, \u+Au' = u" . Multiplying 
by the conjugate u and integrating over R yields 

(7.9) A||uj|| 2 + / ajiUjU^dx + / aj2Uju' 2 dx = / uiu'[dx, j = 1,2, 

' Jul Jr Jm. 

where aji is the j-i entry of A. Integrating the righthand side by parts and rearranging gives 



(7.10) 



1 1 1 2 1 11 / 1 1 2 

\Uj\\ +\\Uj\\ 



ajiUju'idx — / aj2UjU 2 dx, j = 1, 2. 



Taking the real part of (7.10) for K(A) > 0, and applying Cauchy's inequality, we have 
(7.11) 



SR ( A ) 1 1 ii j 1 1 2 + ||u'-|| 2 <C(f \u j \\u' l \dx+ I \uj\\u 2 \dx) < C(2e\\ Uj \\ 2 + ^HH 2 + ^-\\u' 2 \\ 2 ), 



for j = 1,2 where (using sin 28 = 2 cos 9 sin 9) C = max^gjj |aj 5 j(x)| = 1. Similarly, 

1 - ,,,2 1 



(7.12) 



|9(A)|||wj|| < C{2£\\u j \\ I + 



4e' 



u 



+ 



.1 l|2 



4e' 



)■ 



Summing (7.11) and (7.12) for j = 1,2 we have 



Oft(A) + |3(A)|)(||m|| 2 + ||u 2 || 2 ) + IKH 2 + ||4|| 2 < C{Ae\\ Ul \\ 2 + 4e||n 2 || 2 + ^|K|| 2 + ^||n 2 || 2 ). 

Taking e = C, we obtain (K(A) + |3(A)|)(| |«i| | 2 + ||u 2 || 2 ) < 4C 2 (||ui|| 2 + ||u 2 || 2 ), so that for 
K(A) > 0, |A| < 4C 2 < 4. □ 



7.2. Numerical results. Here we describe our numerical studies for the designer system (7.2). 



7.2.1. Comparison of D(0), v, 5. For system (7.7) with A = 0, explicit orthonormal initializing 
bases R± for the Evans function, i.e., bases for the stable (unstable) subspaces of ^,(—00) (-4(oo))> 
are given by 



R- = 



( cos(M77r)\ 
sin(M77r) 
cos(M77r) 

\sin(M77r) J 



(— sin(M77r)\ \ 
cos(M77r) 


V / 



) 
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R+ = 



( sin(M77r) \ / cos(A/77r) \ 

cos(Af77r) — sin(M77r) 

-sin(Af77r) ' 

\-cos(M 7 7r)/ \ J 



Likewise, initializing bases r± for the linearized traveling-wave ODE w' 



Aw of the unstable 
sin(M77r)^ 
cos( 



(stable) subspaces of j4(±oo) are given by r_ = ~ (sin(M77r)), and r + = ~ I ^Qg^jj^-^,^ ) ■ 

compute D(0) and f with these choices, for which ( |5.6| ) becomes essentially a tautology m 
Figure [9] demonstrates the effects of the rotation in solving the Evans function ODE. 



(a) * 






Figure 9. Rotating model with M = 2.7174, 7 = 0.635. (a) Plot of profile = 
— 7tanh(7x/2) against x. (b) Zoomed in plot of components of A(Q(x)) against x. 
(c)-(d) Zoomed in plot of components of the evolved manifolds in the Evans function 
computations for A = 27/2. 




Figure 10. Rotating model, (a) Plot of D(0) (solid black dots) and —5u/4 (red 
circles) against 7 for M = 2.72 demonstrating the identity -D(O) = —v8/A. (b) 
Plot of -D(O) (solid red), 5 (dot-dashed green), and v (dashed blue) against 7 for 
M = 2.72. We mark with an arrow the general region where the Hopf bifurcation 
occurs, (c) Plot of D(Q) (solid red), 5 (dot-dashed green), and v (dashed blue) 
against M with 7 = 0.635. (d) Plot of 8{-y) for M = 2.72. Note that the first root 
of 5 occurs at 5* = I /{AM) ps 0.09. 



7.2.2. Evans function computations. We now detail the results of our numerical study for the 
designer model. For all our Evans function computations we used the integrated form |7.7| of the 
eigenvalue problem. We used MATLABs implicit ODE solver odel5s with relative and absolute 
error bounds set respectively to le-6 and le-8. We used winding number computations and a divide 



and conquer scheme described in |C.2| to find roots of the Evans function. To obtain a global picture 
of roots, we used Lemma 7.3 to ensure that the radius of the contour on which the Evans function 
is computed was sufficiently large to enclose all possible unstable eigenvalues. We also sometimes 
used a smaller radius for speed when looking just for the existence of Hopf bifurcations. 

To search for complex roots and or Hopf bifurcation, we held fixed one of the parameters 7, M, 
and M7 while letting a second parameter vary, examining the resulting motion of eigenvalues along 
these curves. The results of such experiments can be fairly complicated, and can involve a large 
number of roots (up to 16 in the experiment of Figure [TT|b)). However, most of these roots are real 
for most of the time, and, when they do collide and split into a complex pair, tend quickly to rejoin 



Note that vectors (*,*,0,0) T are stationary solutions of \7.7\ , so that the Wronskian of the associated modes 
simplifies to 5 times a 2 x 2 block involving third and fourth components, i.e., the Wronskian v of w' = Aw. 
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and become again real. Thus, to find complex roots, and in particular to identify Hopf bifurcations 
(roots that are not only complex, but cross the imaginary axis as a conjugate pair) involves finding 
a rather small-measure subset of parameter space, and requires a systematic search; they were not 
easy to find! An approach that proved quite useful in identifying Hopf bifurcations was to count 
the number of unstable roots at a selection of mesh points in the 7-M7 plane and look for jumps 
of 2 across boundaries of different regions in the number of roots, then check with a refined study 
that (as generically should hold) these correspond to crossing of the imaginary axis by a non-real 
conjugate pair. 

We note in Figure [To] that for M = 2.72 the first root of 5(7) occurs at <5* = 1/(4M) « 0.09. 
Figure 11 (a) demonstrates that Hopf bifurcations exist and occur roughly periodically in the 
parameter M7, with 7 near special values. We used R = 16 for the contour radius and checked 
that 5(7) had no roots in the vicinity of the jumps by 2 to verify that indeed they correspond 
to Hopf bifurcations. Figure 11 (b) shows that for M7 fixed the number of Hopf bifurcations 
occurring increases as 7 — > 0. Figure [IT] (c) shows a Hopf bifurcation occurring as M increases, for 
7 fixed. Figure 12 (c) shows a Hopf bifurcation occurring as 7 increases, for M fixed. (The latter 
corresponds to variation in shock amplitude for a fixed set of equations, hence is our true interest.) 
Figure 12 (d) shows crossing of the two eigenvalues with largest real part, with a double- multiplicity 
eigenvalue occurring at the point of crossing, violating simplicity of the top eigenvalue. 




Figure 11. Designer system, (a) Color plot of the number of roots present for 
(7, M7 = constant). Arrows indicate where the number of roots jumps by 2 corre- 
sponding to a Hopf bifurcation, (b) Color plot of the number of roots present for 
(7, M7 = constant), (c) Plot in the complex plane of the path two roots of the Evans 
function take as they cross the imaginary axis. Here 7 = 0.65 is fixed and M varies. 
There is a root at about A = 0.09 not in the viewing window. (1) When M = 2.57 
there is one root in the viewing window traveling left, (2) at M = 2.5815 a root 
passes through the origin traveling right , (3) the two roots collide at M ~ 2.585, 
(4) upon collision, the roots split off of the real axis, (5) at M ~ 2.661 the roots 
pass through the imaginary axis. 



Remark 7.4. Vanishing of v appears to be closely related to both Hopf bifurcation and appearance 
of complex roots; hence, unlike the situation for gas dynamics, we have evidently violation of 
transversality of profiles within the parameter regime under study. 



Acknowledgement: The numerical Evans function computations performed in this paper were 
carried out using STABLAB [BHZ2 , a MATLAB-based numerical stability package developed by 
Jeffrey Humpherys together with Barker and Zumbrun. We gratefully acknowledge his contribution. 
Thanks also to Stefano Bianchini, Fedja Nazarov, and Ben Texier for interesting conversations on 
this topic, and in particular to Bianchini for his contribution through system (7.2) |Bi]. The second 
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Figure 12. Designer system. In plots (a)-(c) there is a real root to the right of the 
viewing window, (a) Plot of 7 against root location for M = 3.1. The two small 
modulus roots of the Evans function do not collide and split in the right half plane, 
but presumably they do in the left half-plane, (b) For M = 3.2836 two roots of 
the Evans function collide and split into a non real complex conjugate pair which 
cross the imaginary axis indicating a Hopf bifurcation. The arrows correspond to 
the location of the roots of the Evans function for the following values of 7, (1) 
7 ~ 0.6545, (2) 7 ~ 0.664. (c) For M = 3.51 two roots collide and split into a non 
real complex conjugate pair, but they rejoin in the right half plane. The arrows 
correspond to the location of roots for the following values of 7, (1) 7 « 0.6252, (2) 
7 « 0.629, (3) 7 « 0.6409, (4) 7 = 0.641. (d) Plot of the location of the two largest 
modulus roots of the Evans function against 7 for M = 3.2836. Note that the roots 
pass through each other. 
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Figure 13. For (a)-(c), M m 2.730. For (d)-(f), M » 5.8194, and for (g)-(i), 
M 8.9230. We have in (a), (d), and (g) plots of roots in the complex plane as 7 
varies. In (b), (e), and (h) we plot 7 against location. In (c), (f), and (i) we plot 
7 against location where a closed dot corresponds to the real part of the root and 
an open circle corresponds to the imaginary part. In (i) arrows 1 and 2 emphasize 
respectively where two roots on the real axis collide and split to form a complex 
conjugate pair and then collide again and split along the real axis, similarly as in 
(c) and (f). There are 3 roots in (a)-(c), 4 in (d)-(f), and 5 in (g)-(i). 
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Appendix A. Lopatinski computations for partial equation of state 

To compute the Lopatinski condition in terms of the more general relation p = p(r, e) 
necessarily connected with a full equation of state, rewrite (2.8) for smooth solutions as 

n-v x = 0, 

(A.l) v t +p x = 0, 

e t + pv x = 0, 



not 



so that, for w = (r, v, e) T , A± 



-1 

P 



-c, a 2 = 0, a 3 = c, c = y/—p T + pp e , and 













1, *=l 







\-PJ 









(A.2) 



so that hyperbolicity corresponds to p T — pp e < 0. Computing A\ 



o ^+ 




A + r 2 




4° r 




Using (2.12) to compute [U] = [t](1, —a, —p), we obtain for a 1-shock 



(A.3) 



5 = [r] det 



-a 

-p 



~Pe 

Pt 




where all quantities are evaluated at U + , or, using pp e — p 7 



(A.4) 



6[r] 



(a — c)c 2 + p e avc = (a — c)c 2 +p e c[p], 



we obtain 



recovering (Lop 1 ) (noting by homotopy to the small amplitude case that 5 > ~ stability) in the 
form E| < — — , equivalent to (Lop 1 ) by the relation p e = ff = , where p(r, S) := — e r (r, S). 



es 



This yields the alternate forms % < 2[p] and % > l[p] of (Strong) and (Weak), convenient for 
computations in form (2.8) involving only a pressure law p and not a complete equation of state. 
When e T s < 0, this is equivalent to Majda's condition (2.14) with T = ~ - £ts rewritten as T = rp e . 



Remark A.l. As we have seen, the computations carried out here are simpler, if anything, than 
the computations in entropy coordinates (2.10). However, what is not clear from this formulation 
is why the computations simplify so, involving only the quadratic formula for a 3 x 3 determinant, 
a circumstance that originates from the decoupled form of the equations in entropy variables. 
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Appendix B. Helmholtz energy and associated energy functions 



The completion of a pressure law p = P(T, r) by a compatible energy law e = E(T, r) may be 
obtained more systematically by referring to a related circle of ideas involving the Helmholtz energy 

(B.l) A:=e-TS, 

where e, T, and S have their usual meanings. Here, A is the Legendre transform of e with respect to 
S, T = eg] that is, the relation between e and A is analogous to the relation between Hamiltonian 
and Lagrangian energies in classical mechanics. See |Wlj or |MPj for background discussion. For a 
nice dicscussion of the Legendre transform and relations to the principle of least action, see [W2J . 
The key properties of the Helmoltz energy from our point of view are as follows. 

Lemma B.l. Considering A = A(T,t), we have 

(B.2) p = -A T , S = -A T . 

Proof. By direct computation, At = es(dS/dT) — T(dS/dT) — S = —S, verifying the second 
relation. Similarly, 

A T = e T - T{dS/dr) = e T + es(dS/dr) - T{dS/dr) = e T , 
verifying the first. □ 



Evidently, using (B.2), we may recover A from the pressure law p = p(T, r), up to an arbitrary 
function of T, by antidifferentiation in r. More fundamentally 

(B.3) A = -RTlnQ, 

where Q(T, r) is the (rescaled) partition function coming from the statistical thermodynamic deriva- 
tion of the equation of state. 

Exactness. To check whether a pair of gas laws e = e(T, r) and p = p(T, r) are "exact," that 
is, originate from a complete equation of state A = A(T,t), or, equivalently, e = e(r, S), we may 
note that differentiating S = with respect to r yields S T = hence exactness, S T = pt, is 

equivalent to 

(B.4) Tp T = e T +p. 

We can use this to determine also exactness of a pair of gas laws in the alternative form p = p(r, e) , 
T = T(t, e), using the Implicit Function Theorem to rewrite (B.4) as 

(B.5) Tp e =pf e -f T . 

(Here, we are implicitly assuming T e > 0, so that we may invert T = T(e, r) to solve for e = e(T, r).) 

Hyperbolicity and convexity. The hyperbolicity condition may be deduced using p(S, r) = 
p(r, e s (S, t)) to obtain p T = p T +p T es T = Pt - PtPs = p T - PtPt{9T /dS), or, by S = -A T , hence 
dS/dT = -A TT , 

~ , PtPt A t - A TT A TT 

Pt =Pt + -j" 



Att Att 

Thus, hyperbolicity, p T < 0, holds when Att > and A is convex as a function of (r,T), or else 
Att = —dS/dT = —1/ess < (as usual) and A is not concave . Convexity of e may be computed 
similarly to be equivalent to Att < and es^err — = Atx_ > o, or A TT > and Att < 0. 

Example B.2 (Ideal gas). Here, we have A = —RTlnr — CRT, from which we find as usual 
e(T, t) = |i?T (monatomic gas law, 7 = 5/3.) 
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Example B.3 (van der Waals gas). From the reference [WlJ, we find for T = 2/3, 
(B.6) A = -RT In ((r - b)T 3 ^ - °- 

where the rescaled partition function Q(T,r) satisfies In Q(t,T) = ln(r — b) + + ClnT 3 / 2 , 



C = constant, and S = —At = i?(ln ((r — 6)T 3 / 2 ) + |). Thus, to obtain e, we can just compute 
e(T, t) = A + Ts 



IRT 



Example B.4 (Redlich-Kwong gas). Here |Wlj . A = -RT[\n{j 



ciT- 1 / 2 ln(l + 6/r). 



Appendix C. Numerical Evans function protocol 

For completeness, we briefly describe in this appendix the procedure followed in the numerical 
Evans function studies of part II. For further details, see, e.g., jBHZH [BHRZl IBLeZ| [Z4| IZ5j . 

C.l. Profile and Evans function equations. We first recall the convenient general forms derived 



in |BHLyZl] for profile and Evans function equations of a general class of systems (5.2) arising in 
continuum mechanics: specifically, systems for which B has block-diagonal structure 

(C.l) 5 H=(o b£R r * r , Ro(b)>0, 

and, along the traveling-wave profile w, 

(C.2) det(d/ 1 1 1 — crd/{ ) 1 )(u)) ^ (hyperbolic noncharacteristicity) . 



Example C.l. Navier-Stokes equations ( |6.1[ )-(6.2 ) have form ( |5.2[ ), (C.l) with w = (t,v,T), 

Vfi K 



f°(w) = (r,v,e(r,T)+v 2 /2) T , f\w) = (-u,p(r,T),vp(r,T)) T , and b(w) = r~ l 



Example C.2. Designer system (7.2) has form (5.2), (C.l) with w = (u,v) T = f°(w), r = n 



C.l.l. Traveling-wave system. By the block structure assumption (|C.1|), (|5.3|) decomposes into 

(C.3a) 
(C.3b) 



= fl(w) - fl(U-) - a(fAw) - fi{U-)) 
b(w)w' 2 = fliw) - /*([/_) - a(f°(w) - f$(U-)) 



where f? denotes the fcth block of f 3 . By the Implicit Function Theorem, (C.2) guarantees that 
the first equation may be solved, locally, for w\ as a function of u>2, either analytically (preferable) 
as done here for gas dynamics or numerically, so that the second equation defines a flow in W2- 



Viscous shock profiles are computed using the reduced system (C.3b). 



C.l. 2. Integrated Evans system. Changing to co- moving coordinates x 
steady solution, we obtain the linearized eigenvalue equation 



x 



at in which to is a 



XA°w + {A l w)' = (Bw'Y , 
°— df(w), A 1 w:=df 1 (w)w-adf(w)-&B{w){w,w x ), B := B{w) 



(C.4) 
where 

(C.5) A 

Defining W := A°U, we thus have, integrating (C.4), 
(C.6) XW + A 1 (A 



- x w' '■- 
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B((A V 



or. sot t ing / :— ( ^ ^ )^A°)~ 1 W' ^ € i £" +r - alu = solving for ( J. 0) (.-!") ' IT' in the Ural block of 



(C.6) using again the assumption ( C.2[ ) (see |BHLyZl| ), the first-order integrated Evans syst 



em 
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(C.7) Z' = A int Z, A int :=( -A^^J-i A° 2 - A^A 1 ^ A{ 2 



The integrated Evans function D{\) is computed as a Wronskian of (C.7). 



C.2. Winding number computations. The Evans function is analytic in the closed right-half 
plane where, away from A = 0, zeros of the Evans function match eigenvalues of the system in both 
location and multiplicity. Thus we can determine spectral stability of the system by numerically 
computing the winding number of the Evans function on a semicircle B(0,R) n A) > 0} in the 
right-half plane with R chosen large enough to ensure that any possible unstable eigenvalues are 
inside the semicircle. We determine the value of R either analytically by energy estimates as in 
[BrJ or numerically by convergence to theoretically-predicted high-frequency asymptotics 

(C.8) D(X) ~ C x e c ^ 

as in |HLyZl IBLZl iBLeZ} IBZlj IBZ2j . The former method is carried out on a case- by-case basis, the 



latter method supported automatically in the MATLAB-based numerical Evans function package 
STABLAB |BHZ2j . 



Here, we follow the standard practice of solving, not the eigenvalue equations (C.4), but their 



integrated form (C.6) [Go! IZH} IHuZH IHLZj . thus removing the zero of D occuring at A = as a 



result of translational invariance of the underlying system (5.2) (see (5.4)). That is, we compute 
the interated Evans function D instead of D. This allows us to perform well-conditioned winding 
number computations directly through A = throughout most of the computational domain, 
greatly speeding performance- the exception being at or near the precise parameter values at which 
stability transitions appear, at which this difficulty is unavoidable. 

To locate unstable roots when they occur, we use the method of moments as described in |Bro| 
IBJNRZ[ IBZlj . computing moments of roots within a contour T by the generalized winding number 

computation 2^ rj er int r j = 2^ §r A ^ recovering rj by solution of an resulting polynomial 

equation, and using a divide-and-conquer strategy to keep the number of roots per subcontour 
small enough for optimal accuracy vs. speed (in practice, one or two). 

Standard practice (see, e.g., [BHRZl IHLZl IBLeZj ). is to approximate the traveling wave profile 
using MATLAB's boundary- value solver bvp6c [HMj, an adaptive Lobatto quadrature scheme, 
with projective boundary conditions, and on a truncated, finite computational domain [— L, L], 
where the values of approximate plus and minus spatial infinity L are determined experimentally by 
the requirement that the absolute error between the computed profile values at x = ±L be within a 
prescribed tolerance of the actual endstates U± at plus and minus spatial infinity. See |Bel|IBe2] for 
theoretical justification and rigorous error bounds for this approach. We follow this approach for the 



designer system (7.2). However, for our gas-dynamical computations near the large-amplitude limit 
5- — > —00, this becomes impractical due to stiffness/presence of multiple scales in the traveling- 
wave ODE, and we depart from this method in favor of a simple shooting algorithm with a stiff 
ODE solver (specifically, the adaptive mesh solver odel5s supported in MATLAB). 

We then approximate the Evans function D numerically following the approach described in 
detail in [B HZll IBLeZ} IBLZj . To preserve analyticity of the Evans function, we use the method 
of Kato |BHZU IBrZ} IGZl IHuZlt |Kj IZ4} IZ7| to obtain a holomorphic initializing basis with respect 
to A. To evolve in x the manifolds whose determinant yields the Evans function, we use the polar 
coordinate ("analytic orthogonalization" ) method described in [HuZl]. This method evolves an 

41 



orthonormal solution basis along with a complex radial equation that maintains the property of 
analyticity. See [H uZl| IZ7| for theoretical justification and error bounds. To ensure an accurate 
winding number count, an adaptive mesh in A is used requiring that the relative error of change 
in D for each step be less than 0.2. Recall, by Rouche's Theorem, that if the relative variation of 
D is less than 1.0, then winding number accuracy is preserved. These routines are supported in 
STABLAB |BHZ2| . 

For gas-dynamical systems (4.1) and (4.3), the Evans function is poorly conditioned due to 
multiple scales/stiffness of the eigenvalue ODE. In this case, we find it necessary to use MATLAB's 
adaptive mesh stiff differential equation solver odel5s both to solve for the profile and to perform 
the Evans function computations, rather than the Runge-Kutta-Fehlberg method encoded in ode45 
that is typically used in Evans function computations. 



C.3. Hardware and computational statistics. In this short section we detail computational 
statistics for some typical values in the numerical studies reported in previous sections. All numer- 
ical computations were performed on either a Mac Pro with 2 Quad-Core Intel Xeon processors 
with speed 2.26 GHz, the super computer Quarry at Indiana University, or a MacBook with an 
Intel Core 2 Duo processor with speed 2.0 GHz. 

The Evans function computations were carried out using MATLAB's implicit ODE solver, odel5s, 
with relative error tolerance set at le-6 and absolute error set at le-8. For most of the parameters 
in our studies, we used MATLAB's boundary value solver bvp5c with relative error set at le-6, 
absolute error set at le-8, and the boundary error set at le-6. For the few parameters for which 
the boundary value solver failed, we used a shooting method as described above, with odel5s. 

We now relate computational statistics, run on the MacBook, for each example system. In each 
case, a minimum of 40 A points were used on a semi-circular mesh with radius R = 10. Additional 
points were added to the mesh adaptively so that the relative change in consecutive contour points 
varied by no more than 0.2. The full contour was constructed by using the conjugate symmetry 
of the Evans function. For C = 10 and S- = —4 in the global model, it took 23.6 seconds to 
compute the profile using a boundary value solver with a tanh solution as an initial guess. It took 
176 seconds compute the Evans function on 52 points. For S- = —4 in the local model it took 
62.1 seconds to compute the profile using a boundary value solver and a tanh solution for an initial 
guess. It took 95.5 seconds to compute the Evans function on 63 points. For S- = —4 in the stable 
model it took 7.98 seconds to solve the profile with a boundary value solver initialized with a tanh 
solution. It took 28.9 seconds to compute the Evans function on 40 points. 

To give some perspective on the computational difficulty of the example models, we include data 
for the isentropic gas model considered in [BHRZJ with 7 = 5/3 and v + = le — 4. The boundary 
value solver took 0.899 seconds using a tanh solution for an initial guess. The Evans function 
took 27.7 seconds to compute on 40 points. As the above data indicates, solving the profile in the 
example gas systems considered here is much more difficult than in isentropic gas. The profile was 
amenable to solve using a boundary value solver for S- as small as S- = —200 for the global model 
with C = 10, S- = —8.1 for the local model, and S- = —11.5 for the stable model. When using 
continuation to solve, we took steps of 0.1 in S-. 

We also compare computational performance on the Mac Pro with its 8 cores, which is able to 
utilize parallel computing. For C = 10 and S- = —4 in the global model, it took 16.9 seconds 
to compute the profile using a boundary value solver with a tanh solution as an initial guess. It 
took 28.75 seconds compute the Evans function on 52 points. For S- = —4 in the local model it 
took 37.78 seconds to compute the profile using a boundary value solver and a tanh solution for an 
initial guess. It took 17.97 seconds to compute the Evans function on 63 points. For S- = —4 in 
the stable model it took 5.63 seconds to solve the profile with a boundary value solver initialized 
with a tanh solution. It took 6.81 seconds to compute the Evans function on 40 points. For the 
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isentropic model with 7 = 5/3 and v + = le — 4, the boundary value solver took 0.5833 seconds 
using a tanh solution for an initial guess. The Evans function took 6.0929 seconds to compute on 
40 points. 
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